Topology Optimization for FDM Parts Considering the Hybrid Deposition Path Pattern

Abstract

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.

1. 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].

Topology optimization has been widely treated as the main computational design method for AM [1,2,3,4], because it explores a large design space and successfully applies in different physical disciplines [5,6,7,8]. Diverse topology optimization methods have been proposed, including the homogenization method [9,10], solid isotropic material with penalization (SIMP) [11,12], evolutionary structural optimization (ESO) [13], the level set method (LSM) [14,15], the moving morphable component (MMC) method [16], and some others.

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 profile offsets through a constant distance $T$, where the deposition path is assumed to be distributing along the outer profile and the interior substrate structure is filled with fixed deposition paths of angle $\overline{\theta }$ (Figure 1). Under this scenario, to better close the practice, a level set-based method was proposed by Liu et al. to take the HDP pattern into account to address the material anisotropy topology optimization [28]. Besides, the field of fiber reinforcement composites is also highly relevant to additive manufacturing-oriented topology optimization [29,30]. An interdependent two-level optimization approach was proposed by Humberto et al. to optimize both fiber angle and intrinsic thickness, and producible results could be obtained by this method [31]. A novel framework for the optimized topology and the fiber paths was developed in [32]; the optimized distribution of the material and the fiber orientation are achieved by two methods: a density-based method and the level set method.

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.

2. 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.

2.1. The Formulation for the Boundary Layer

According to the classical laminate theory, ${\mathit{D}}_0$ is the laminae unrotated compliance tensor [39,40]:

$${\mathit{D}}_0=\left[\begin{array}{ccc}\frac{E_x}{1-v_{xy}{v}_{yx}}& \frac{v_{yx}{E}_x}{1-v_{xy}{v}_{yx}}& 0\\ \frac{v_{yx}{E}_x}{1-v_{xy}{v}_{yx}}& \frac{E_y}{1-v_{xy}{v}_{yx}}& 0\\ 0& 0& G_{xy}\end{array}\right],$$

(1)

and $\mathit{D}\left(\theta \right)$ is the 2D orthotropic elasticity tensor given any angle $\theta$:

$$\mathit{D}\left(\theta \right)=\mathit{T}\left(\theta \right){\mathit{D}}_0\mathit{T}{\left(\theta \right)}^T,$$

(2)

with $\mathit{T}\left(\theta \right)$ being the transform matrix which is used to conduct the matrix coordinate transform:

$$\mathit{T}\left(\theta \right)=\left[\begin{array}{ccc}{\cos }^2\theta & {\sin }^2\theta & -2\sin \theta \cos \theta \\ {\sin }^2\theta & {\cos }^2\theta & 2\sin \theta \cos \theta \\ \sin \theta \cos \theta & -\sin \theta \cos \theta & {\cos }^2\theta -{\sin }^2\theta \end{array}\right]$$

(3)

In this work, the local fiber orientation $\theta$ could be analytically expressed by Equation (3) and counted in the counter-clockwise direction, as shown in Figure 3.

$$\mathit{\theta }=\frac{\pi }{2}+arctan\left(\frac{\frac{\partial \widehat{\mathit{\phi }}}{\partial y}}{\frac{\partial \widehat{\mathit{\phi }}}{\partial x}}\right)$$

(4)

$\frac{\partial \widehat{\mathit{\phi }}}{\partial y}$ and $\frac{\partial \widehat{\mathit{\phi }}}{\partial x}$ are the spatial gradients of the filtered design field $\widehat{\mathit{\phi }}$, and their derivations could refer to later content.

Thus, the element stiffness matrix could be written as:

$$\mathit{K}\mathit{e}=\delta \left(\mathit{D}\right)$$

(5)

where $\delta (*)$ is the element stiffness matrix assembly operator.

2.2. The Optimization Model

In this paper, a standard compliance minimization problem, subject to mass constraint, is studied. The optimization problem is formulated as:

$$\{\begin{array}{c}Find:{\mu }_E\left(E=1,2,\dots ,N\right)\\ \begin{array}{c}Minimize:C={\mathit{U}}^T\mathit{K}\mathit{U}\end{array}\\ Subjectto:\{\begin{array}{c}\mathit{K}\mathit{U}=\mathit{F}\\ G\le M_d\\ 0<{\mu }_{min}\le {\mu }_E\le 1\end{array}\end{array}$$

(6)

where ${\mu }_E$ is the design variable and the Eth element density and $N$ is the total number of elements. $\mathit{K}$, $\mathit{U}$, and $\mathit{F}$ are the global stiffness matrix, displacement vector, and force vector, respectively. ${\mu }_{min}$ is a small number to avoid matrix singularity. $G$ is the mass fraction constrained by the maximum mass fraction $M_d$. Note that the detail expression of $G$ could refer to later content.

2.3. Material Interpolation Strategy

Based on the SOMP, the element stiffness matrix and physical density interposition of the element $E$ are defined as an interpolation of ${\phi }_E$ and $\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}$:

$$\begin{array}{c}\begin{array}{c}{\mathit{K}}_E\left({\phi }_E,\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)={\left({\phi }_E\right)}^p·\mathit{K}{\mathit{e}}_2+{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^p·\overline{\mathit{K}{\mathit{e}}_{1,E}}-{\left({\phi }_E\right)}^p·{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^p·\mathit{K}{\mathit{e}}_2\\ {\rho }_E\left({\phi }_E,\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)={\rho }_S·{\phi }_E+{\rho }_B·\left(1-{\rho }_S·{\phi }_E\right)·{\left(\overline{{\Vert \nabla \widehat{\phi }\Vert }_{\alpha }}\right)}_E\end{array}\end{array}$$

(7)

where the $\overline{\mathit{K}{\mathit{e}}_{1,E}}$ is the stiffness matrix of the mth element that belongs to the boundary layer, $\mathit{K}{\mathit{e}}_2$ is the element stiffness matrix for the substrate material. $p$ is employed to penalize the intermediate densities, so as to derive the black and white solution. Note that when $\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}$ approaches zero, i.e., when going away from the boundary layer, the expressions reduce to:

$$\begin{array}{c}\begin{array}{c}{\mathit{K}}_E\left({\phi }_E,0\right)={\left({\phi }_E\right)}^p·\mathit{K}{\mathit{e}}_2\\ {\rho }_E\left({\phi }_E,0\right)={\rho }_S·{\phi }_E\end{array}\end{array}$$

(8)

where ${\rho }_S$ is the element density in the substrate domain.

Meanwhile, at the other extreme, when $\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}$ approaches $1$, i.e., at the boundary layer region, the expressions could be expressed as:

$$\begin{array}{c}\begin{array}{c}{\mathit{K}}_E\left({\phi }_E,1\right)=\overline{\mathit{K}{\mathit{e}}_{1,E}}\\ {\rho }_E\left({\phi }_E,1\right)={\rho }_B\end{array}\end{array}$$

(9)

where ${\rho }_B$ is the element density in the boundary area.

In conclusion, both the physical density and stiffness are interpolated based on the projected design variable field, ${\phi }_E$, and the normalized gradient of the filtered field $\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}$.

2.4. Objective Function

The structural compliance is equal to the sum of the element strain energy, which could be expressed as:

$$C={\displaystyle \sum }_{E=1}^N{Q}_E$$

(10)

where $Q_E$ is the strain energy of the element $E$:

$$Q_E={\mathit{U}}_E{}^T{\mathit{K}}_E{\mathit{U}}_E$$

(11)

Substituting the element stiffness matrix interpolation in Equation (7) into Equation (11) will yield:

$$\begin{array}{c}{Q}_E={\left({\phi }_E\right)}^p·\left({\mathit{U}}_E{}^T\mathit{K}{\mathit{e}}_2{\mathit{U}}_E\right)+{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^p·\left({\mathit{U}}_E{}^T\overline{\mathit{K}{\mathit{e}}_{1,E}}{\mathit{U}}_E\right)\\ -{\left({\phi }_E\right)}^p·{\left(\overline{{\Vert \nabla \widehat{\phi }\Vert }_{\alpha }}\right)}_E{}^p·\left({\mathit{U}}_E{}^T\mathit{K}{\mathit{e}}_2{\mathit{U}}_E\right)\end{array}$$

(12)

For brevity, the following notation is introduced:

$$Q_E={\omega }_S+{\omega }_G+{\omega }_{SG}={\epsilon }_E^S·Q{2}_E+{\epsilon }_E^G·Q{1}_E-{\epsilon }_E^{SG}·Q{2}_E$$

(13)

where ${\epsilon }_E^S$, ${\epsilon }_E^G$, and ${\epsilon }_E^{SG}$ could be treated as three penalized pseudo-density fields:

$$\{\begin{array}{c}{\epsilon }_E^S={\left({\phi }_E\right)}^p\\ {\epsilon }_E^G={\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^p\\ {\epsilon }_E^{SG}={\left({\phi }_E\right)}^p·{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^p\end{array}$$

(14)

$Q{2}_E$ is the element strain energy for the substrate element and $Q{1}_E$ is treated as the modified element strain energy for the boundary layer element:

$$\{\begin{array}{c}Q{1}_E=U_E{}^T\overline{\mathit{K}{\mathit{e}}_{1,E}}{U}_E\\ Q{2}_E=U_E{}^T\mathit{K}{\mathit{e}}_2{U}_E\end{array}$$

(15)

2.5. Mass Constraint

In this paper, the mass constraint function could be written as:

$$G={\displaystyle \sum }_{E=1}^N{G}_E$$

(16)

where $G_E$ is the mass of the element $E$:

$$G_E=M_0·{\rho }_E=M_0·\left[{\rho }_S·{\phi }_E+{\rho }_B·\left(1-{\rho }_S·{\phi }_E\right)·{\left(\overline{{\Vert \nabla \widehat{\phi }\Vert }_{\alpha }}\right)}_E\right]$$

(17)

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:

$$G_E=M_0·\left({ϵ}_E^S+{ϵ}_E^G-{ϵ}_E^{SG}\right)$$

(18)

where $M_0$ is the design element standard mass; ${ϵ}_E^S$, ${ϵ}_E^G$, and ${ϵ}_E^{SG}$ are the three non-penalized pseudo-design fields, as follows:

$$\{\begin{array}{c}{ϵ}_E^S={\rho }_S·{\phi }_E\\ {ϵ}_E^G={\rho }_B·\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\\ {ϵ}_E^{SG}={\rho }_S·{\rho }_B·{\phi }_E·\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\end{array}$$

(19)

2.6. 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.

2.6.1. Sensitivity Analysis for Objective Function

Recalling Equation (13), $\frac{\partial Q_E}{\partial {\mu }_E}$ could be written in the following form:

$$\frac{\partial Q_E}{\partial {\mu }_E}=\frac{\partial {\omega }_S}{\partial {\mu }_E}+\frac{\partial {\omega }_G}{\partial {\mu }_E}+\frac{\partial {\omega }_{SG}}{\partial {\mu }_E}$$

(20)

The first derivative term could be obtained using the chain rule:

$$\frac{\partial {\omega }_S}{\partial {\mu }_E}=\frac{\partial {\epsilon }_E^S}{\partial {\mu }_E}·Q{2}_E=\frac{\partial {\left({\phi }_E\right)}^p}{\partial {\mu }_E}·Q{2}_E=p·{\phi }_E{}^{p-1}·\frac{\partial {\phi }_E}{\partial {\widehat{\mu }}_E}·\frac{\partial \widehat{{\mu }_E}}{\partial {\mu }_E}·Q{2}_E$$

(21)

The term $\frac{\partial {\phi }_E}{\partial {\widehat{\mu }}_E}$ represents the standard modification of sensitivities due to projection and the detail expression can be found in [41]. $\frac{\partial \widehat{{\mu }_E}}{\partial {\mu }_E}$ 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:

$$\frac{\partial {\omega }_G}{\partial {\mu }_E}=\frac{\partial {\epsilon }_E^G}{\partial {\mu }_E}·Q{1}_E+\frac{\partial Q{1}_E}{\partial {\mu }_E}·{\epsilon }_E^G$$

(22)

and the first derivative term could be written as the following form:

$$\frac{\partial {\epsilon }_E^G}{\partial {\mu }_E}=\frac{\partial {\epsilon }_E^G}{\partial {\widehat{\phi }}_N}·\frac{\partial {\widehat{\phi }}_N}{\partial {\phi }_N}·\frac{\partial {\phi }_N}{\partial {\phi }_E}·\frac{\partial {\phi }_E}{\partial {\widehat{\mu }}_E}·\frac{\partial {\widehat{\mu }}_E}{\partial {\mu }_E}$$

(23)

where

$$\frac{\partial {\epsilon }_E^G}{\partial {\widehat{\phi }}_N}=p·{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^{p-1}·\frac{\partial \left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}{\partial {\widehat{\phi }}_N}=p·{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^{p-1}·\frac{\partial \left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}·\frac{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}{\partial {\widehat{\phi }}_N}$$

(24)

in which the term $\frac{\partial \left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}$ indicates the standard modification of sensitivities due to the projection of the ${\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }$ field. The derivative of the normalized gradient norm $\frac{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}{\partial {\widehat{\phi }}_N}$ could be written as:

$$\frac{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}{\partial {\widehat{\phi }}_N}=\frac{\alpha }{\Vert \nabla \widehat{\phi }\Vert }·\left(\frac{\partial {\widehat{\phi }}_E}{\partial x}\frac{\partial \frac{\partial {\widehat{\phi }}_E}{\partial x}}{\partial {\widehat{\phi }}_N}+\frac{\partial {\widehat{\phi }}_E}{\partial y}\frac{\partial \frac{\partial {\widehat{\phi }}_E}{\partial y}}{\partial {\widehat{\phi }}_N}\right)$$

(25)

in which:

$$\begin{array}{c}\frac{\partial {\widehat{\phi }}_E}{\partial x}=\frac{\partial \left(N^T{\widehat{\phi }}_N\right)}{\partial x}=B_x{\widehat{\phi }}_N\\ \frac{\partial {\widehat{\phi }}_E}{\partial y}=\frac{\partial \left(N^T{\widehat{\phi }}_N\right)}{\partial y}=B_y{\widehat{\phi }}_N\end{array}$$

(26)

where $\mathit{N}$ is a vector of the four shape functions relating nodal variable ${\widehat{\phi }}_N$ with the elemental variable ${\widehat{\phi }}_E$; $B_x$ and $B_y$ are the gradient computation matrices for $\mathit{N}$ in the $x$ and $y$ directions, respectively, and are independent from the design variables. Therefore, we have:

$$\begin{array}{c}\frac{\partial \frac{\partial {\widehat{\phi }}_E}{\partial x}}{\partial {\widehat{\phi }}_N}=\frac{\partial \left(B_x{\widehat{\phi }}_N\right)}{\partial {\widehat{\phi }}_N}=B_x\\ \frac{\partial \frac{\partial {\widehat{\phi }}_E}{\partial y}}{\partial {\widehat{\phi }}_N}=\frac{\partial \left(B_y{\widehat{\phi }}_N\right)}{\partial {\widehat{\phi }}_N}=B_y\end{array}$$

(27)

Thus,

$$\frac{\partial {\epsilon }_E^G}{\partial {\widehat{\phi }}_N}=p·{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^{p-1}·\frac{\partial \left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}·\frac{\alpha }{\Vert \nabla \widehat{\phi }\Vert }·\left(\frac{\partial {\widehat{\phi }}_E}{\partial x}{\mathit{B}}_x+\frac{\partial {\widehat{\phi }}_E}{\partial y}{\mathit{B}}_y\right)$$

(28)

and the following notation of Equation (28) is introduced:

$$\frac{\partial {\epsilon }_E^G}{\partial {\widehat{\phi }}_N}=M_E{}^x{\mathit{B}}_x+M_E{}^y{\mathit{B}}_y$$

(29)

where

$$\begin{array}{c}{M}_E{}^x=p·{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^{p-1}·\frac{\partial \left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}·\frac{\alpha }{\Vert \nabla \widehat{\phi }\Vert }·\frac{\partial {\widehat{\phi }}_E}{\partial x}\\ M_E{}^y=p·{\left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}^{p-1}·\frac{\partial \left(\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\right)}{\partial \left({\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }\right)}·\frac{\alpha }{\Vert \nabla \widehat{\phi }\Vert }·\frac{\partial {\widehat{\phi }}_E}{\partial y}\end{array}$$

(30)

For the term $\frac{\partial Q{1}_E}{\partial {\mu }_E}$, it could have following form:

$$\frac{\partial Q{1}_E}{\partial {\mu }_E}=\frac{\partial Q{1}_E}{\partial {\theta }_E}·\frac{\partial {\theta }_E}{\partial {\widehat{\phi }}_N}·\frac{\partial {\widehat{\phi }}_N}{\partial {\phi }_N}·\frac{\partial {\phi }_N}{\partial {\phi }_E}·\frac{\partial {\phi }_E}{\partial {\widehat{\mu }}_E}·\frac{\partial {\widehat{\mu }}_E}{\partial {\mu }_E}$$

(31)

where

$$\frac{\partial Q{1}_E}{\partial {\theta }_E}=\frac{\partial \left({\mathit{U}}_E{}^T\mathit{K}{\mathit{e}}_{1,E}{\mathit{U}}_E\right)}{\partial {\theta }_E}={\mathit{U}}_E{}^T\delta \left(\frac{\partial {\mathit{D}}_E}{\partial {\theta }_E}\right){\mathit{U}}_E$$

(32)

Substituting the material interpolation of Equation (2) into the term $\frac{\partial {\mathit{D}}_E}{\partial {\theta }_E}$ will yield:

$$\frac{\partial {\mathit{D}}_E}{\partial {\theta }_E}={\mathit{U}}_E{}^T\delta \left(\frac{\partial \left(\mathit{T}\left(\theta \right){\mathit{D}}_0\mathit{T}{\left(\theta \right)}^T\right)}{\partial {\theta }_E}\right){\mathit{U}}_E$$

(33)

The elastic tensor of matrix material ${\mathit{D}}_0$ is independent of ${\mu }_E$. Using the chain rule again, we could arrive at:

$$\frac{\partial \left(\mathit{T}\left(\theta \right){\mathit{D}}_0\mathit{T}{\left(\theta \right)}^T\right)}{\partial {\theta }_E}=\frac{\partial \mathit{T}\left({\theta }_E\right)}{\partial {\theta }_E}{\mathit{D}}_0\mathit{T}{\left({\theta }_E\right)}^T+\mathit{T}\left({\theta }_E\right){\mathit{D}}_0\frac{\partial \mathit{T}{\left({\theta }_E\right)}^T}{\partial {\theta }_E}$$

(34)

where

$$\frac{\partial \mathit{T}\left({\theta }_E\right)}{\partial {\theta }_E}=\left[\begin{array}{ccc}-2·\sin \theta ·\cos \theta & \sin 2\theta & -2·\cos 2\theta \\ \sin 2\theta & -2\sin \theta ·\cos \theta & 2·\cos 2\theta \\ \cos 2\theta & -\cos 2\theta & -2·\sin \theta ·\cos \theta -\sin 2\theta \end{array}\right]$$

(35)

Then, the derivative of $\frac{\partial {\theta }_E}{\partial {\widehat{\phi }}_N}$ is obtained as:

$$\frac{\partial {\theta }_E}{\partial {\widehat{\phi }}_N}=\frac{1}{1+{\left(\frac{\partial {\widehat{\phi }}_E}{\partial y}\right)}^2·{\left(\frac{\partial {\widehat{\phi }}_E}{\partial x}\right)}^{-2}}·\frac{\frac{\partial {\widehat{\phi }}_E}{\partial x}·{\mathit{B}}_y-\frac{\partial {\widehat{\phi }}_E}{\partial y}·{\mathit{B}}_x}{{\left(\frac{\partial {\widehat{\phi }}_E}{\partial x}\right)}^2}$$

(36)

Similarly, the following notation of Equation (36) is introduced:

$$\frac{\partial {\theta }_E}{\partial {\widehat{\phi }}_N}={{\overline{M}}_E}^x{\mathit{B}}_y+{{\overline{M}}_E}^y{\mathit{B}}_x$$

(37)

where

$$\begin{array}{c}{{\overline{M}}_E}^x=\frac{1}{1+{\left(\frac{\partial {\widehat{\phi }}_E}{\partial y}\right)}^2·{\left(\frac{\partial {\widehat{\phi }}_E}{\partial x}\right)}^{-2}}·\frac{1}{\frac{\partial {\widehat{\phi }}_E}{\partial x}}\\ {{\overline{M}}_E}^y=\frac{1}{1+{\left(\frac{\partial {\widehat{\phi }}_E}{\partial y}\right)}^2·{\left(\frac{\partial {\widehat{\phi }}_E}{\partial x}\right)}^{-2}}·\frac{\frac{\partial {\widehat{\phi }}_E}{\partial y}}{{\left(\frac{\partial {\widehat{\phi }}_E}{\partial x}\right)}^2}\end{array}$$

(38)

The third derivative term of Equation (20) could be obtained simply by using the product rule and all the terms already obtained above:

$$\frac{\partial {\epsilon }_E^{SG}}{\partial {\mu }_E}·Q{2}_E=\frac{\partial {\epsilon }_E^S}{\partial {\mu }_E}·{\epsilon }_E^G·Q{2}_E+\frac{\partial {\epsilon }_E^G}{\partial {\mu }_E}·{\epsilon }_E^S·Q{2}_E$$

(39)

2.6.2. Sensitivity Analysis for Mass Constraint

For the term $\frac{\partial G_E}{\partial {\mu }_E}$, it could be similarly elaborated by the following expressions:

$$\frac{\partial G_E}{\partial {\mu }_E}=M_0·\frac{\partial {\rho }_E}{\partial {\mu }_E}=M_0·\frac{\partial \left({ϵ}_E^S+{ϵ}_E^G-{ϵ}_E^{SG}\right)}{\partial {\mu }_E}$$

(40)

where,

$$\{\begin{array}{c}\frac{\partial {ϵ}_E^S}{\partial {\mu }_E}={\rho }_S·\frac{\partial {\phi }_E}{\partial \widehat{{\mu }_E}}·\frac{\partial \widehat{{\mu }_E}}{\partial {\mu }_E}\\ \frac{\partial {ϵ}_E^G}{\partial {\mu }_E}={\rho }_B·\frac{\partial \overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}}{\partial {\mu }_E}\\ \frac{\partial {ϵ}_E^{SG}}{\partial {\mu }_E}={\rho }_S·{\rho }_B·\left(\frac{\partial {ϵ}_E^S}{\partial {\mu }_E}·{ϵ}_E^G+\frac{\partial {ϵ}_E^G}{\partial {\mu }_E}·{ϵ}_E^S\right)\end{array}$$

(41)

All the derivative terms in Equation (41) could be obtained by the previous expressions.

2.6.3. 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]:

$$-R^2{\nabla }^2\widehat{\mathit{\mu }}+\widehat{\mathit{\mu }}=\mathit{\mu }$$

(42)

by imposing the homogeneous Neumann boundary conditions ($\frac{\partial \widehat{\mathit{\mu }}}{\partial \mathit{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), $\mathit{\mu }$ represents the unfiltered design field, and $\widehat{\mathit{\mu }}$ 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]:

$$R=\frac{r_{min}}{2\sqrt{3}}$$

(43)

Using the finite element method to discrete Equation (42):

$${\mathit{K}}_x{\widehat{\mathit{\mu }}}_N={\mathit{\mu }}_N$$

(44)

where ${\mathit{K}}_x$ is the standard stiffness matrix in finite element method for the scalar problem corresponding to $R_x$ and ${\widehat{\mathit{\mu }}}_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:

$$\frac{\partial {\widehat{\mathit{\mu }}}_N}{\partial {\mathit{\mu }}_N}={\mathit{K}}_x{}^{-1}$$

(45)

The element node representation of the field is obtained by:

$${\mathit{\mu }}_N={\mathit{T}}_F\mathit{\mu }$$

(46)

where ${\mathit{T}}_F$ is a matrix which maps the elemental values $\mathit{\mu }$ to a vector with nodal values ${\mathit{\mu }}_N$. Similarly, the element node representation of the filtered field could be expressed as:

$$\widehat{\mathit{\mu }}={\mathit{T}}_F{}^T{\widehat{\mathit{\mu }}}_N$$

(47)

and

$$\frac{\partial {\mathit{\mu }}_N}{\partial \mathit{\mu }}={\mathit{T}}_F$$

(48)

According to the above derivation, the filtered sensitivities of the term ${\mathit{\omega }}_S$ can be calculated as:

$$\frac{\partial {\mathit{\omega }}_S}{\partial \mathit{\mu }}={\mathit{T}}_F{}^T\left({\mathit{K}}_1{}^{-1}\left({\mathit{T}}_F\left(p·{\mathit{\phi }}^{p-1}·\frac{\partial \mathit{\phi }}{\partial \widehat{\mathit{\mu }}}·\mathit{Q}\mathbf{2}\right)\right)\right)$$

(49)

the filtered sensitivities of the term ${\mathit{\omega }}_G$ could be rewritten as:

$$\frac{\partial {\mathit{\omega }}_G}{\partial \mathit{\mu }}=\frac{\partial {\mathit{\omega }}_{G1}}{\partial \mathit{\mu }}+\frac{\partial {\mathit{\omega }}_{G2}}{\partial \mathit{\mu }}$$

(50)

The terms $\frac{\partial {\mathit{\omega }}_{G1}}{\partial \mathit{\mu }}$ and $\frac{\partial {\mathit{\omega }}_{G2}}{\partial \mathit{\mu }}$ can be computed as:

$$\begin{array}{c}\frac{\partial {\mathit{\omega }}_{G1}}{\partial \mathit{\mu }}={\mathit{T}}_F{}^T\left({\mathit{K}}_1{}^{-1}\left({\mathit{T}}_F\left(\mathit{d}\mathit{c}·\frac{\partial \mathit{\phi }}{\partial \widehat{\mathit{\mu }}}\right)\right)\right)\\ \mathit{d}\mathit{c}={\mathit{T}}_F{}^T\left({\mathit{K}}_2{}^{-1}\left({\mathit{T}}_F{\left({\mathit{M}}^x·\mathit{Q}\mathbf{1}\right)}^T{\mathit{B}}_X{}^T+{\left({\mathit{M}}^y·\mathit{Q}\mathbf{1}\right)}^T{\mathit{B}}_Y{}^T\right)\right)\end{array}$$

(51)

and

$$\begin{array}{c}\frac{\partial {\mathit{\omega }}_{G2}}{\partial \mathit{\mu }}={\mathit{T}}_F{}^T\left({\mathit{K}}_1{}^{-1}\left({\mathit{T}}_F\left(\overline{\mathit{d}\mathit{c}}·\frac{\partial \mathit{\phi }}{\partial \widehat{\mathit{\mu }}}\right)\right)\right)\\ \overline{\mathit{d}\mathit{c}}={\mathit{T}}_F{}^T\left({\mathit{K}}_2{}^{-1}\left({\mathit{T}}_F\left({\left({\overline{\mathit{M}}}^x·{\mathit{\epsilon }}^G\right)}^T{\mathit{B}}_X{}^T+{\left({\overline{\mathit{M}}}^y·{\mathit{\epsilon }}^G\right)}^T{\mathit{B}}_Y{}^T\right)\right)\right)\end{array}$$

(52)

Again, the filtered derivative for the term ${\mathit{\omega }}_{SG}$ and the filtered sensitivity analysis for the mass constraint could be obtained in a similar way and are therefore omitted here.

3. 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 ${\beta }_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 ${\beta }_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 ${\mu }_S$ and ${\mu }_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.

4. Case Studies

4.1. Messerschmidt–Bölkow–Blohm (MBB) Problem

4.1.1. The Fully Infilled Substate Problem

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.

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 (${\rho }_S={\rho }_B$) and the direction and the rotation angle $\theta$ 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, ${\rho }_E\ge 0.95$ indicates a clearly formed substrate domain, and $\overline{{\Vert \nabla \widehat{{\phi }_E}\Vert }_{\alpha }}\ge 0.5$ represents a clearly formed boundary layer. The final compliance is $88.3850$, and the optimization terminates at the $320$th 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 ${\beta }_S$ and ${\beta }_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 $200$th 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].

4.1.2. 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, ${\rho }_E\ge 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.

4.1.3. 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 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.

4.1.4. The Influence of Different Boundary Layer Widths

In Figure 11, the same design problem as above is solved with the raster direction $\theta =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.

4.1.5. The Mesh Independence

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.

4.2. Cantilever Problem

4.2.1. 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\times 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.

4.2.2. 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\times 50$) adopted in the same case in Figure 15, it is impossible to identify a boundary layer that is equal or above $T=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.

4.3. Short Cantilever Problem

4.3.1. 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\times 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.

4.3.2. 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, and the derived objective values are $41.4360$, $14.6845$, 27.5177, and 28.1572 for wiggle, honeycomb, $45°$ line, and $-45°$ line, respectively.

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\times 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 2$00\times 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.

5. 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.