Topology Optimization for Multipatch Fused Deposition Modeling 3D Printing

Featured Application: This work can be used for lightweight design of drone structures. Abstract: This paper presents a hybrid topology optimization method for multipatch fused deposition modeling (FDM) 3D printing to address the process-induced material anisotropy. The ‘multipatch’ concept consists of each printing layer disintegrated into multiple patches with di ﬀ erent zigzag-type ﬁlament deposition directions. The level set method was employed to represent and track the layer shape evolution; discrete material optimization (DMO) model was adopted to realize the material property interpolation among the patches. With this set-up, a concurrent optimization problem was formulated to simultaneously optimize the topological structure of the printing layer, the multipatch distribution, and the corresponding deposition directions. An asynchronous starting strategy is proposed to prevent the local minimum solutions caused by the concurrent optimization scheme. Several numerical examples were investigated to verify the e ﬀ ectiveness of the proposed method, while satisfactory optimization results have been derived.


Introduction
In additive manufacturing (AM) or 3D printing, many of the design complexity constraints in conventional manufacturing methods can be eliminated due to its layer-by-layer material deposition nature. Consequently, the creativity of the designers can be greatly released and realized. This is why AM has attracted so much attention in research and industrial applications in the past decade. At present, intense research on AM is being carried out regarding geometry design, material design, computation tools, and manufacturing tools and process development [1,2].
Topology optimization has been treated as the main computational design method for AM [1][2][3]. Because the shape and topology are concurrently designed, topology optimization makes the greatest design freedom possible compared with shape-only or size-only optimization. However, there are new design rules and unique constraints induced by AM, which introduce new challenges such as support structure design/elimination [4][5][6][7][8][9][10], minimum component size constraints [11][12][13][14][15][16], directional material properties [17][18][19], topology design interpretation [20][21][22][23], variable-density cellular structure design [24][25][26][27][28], and many others [29,30]. Detailed reviews on these topics can be found in [2,31]. Despite the efforts, a lack of solutions is still the common issue in the topology optimization for AM as summarized by the literature. Among these topics, the design constraints imposed by the anisotropic properties of AM materials are of interest and will be addressed by the hybrid topology optimization technique proposed in this work.
The directional material properties are rooted in the layer-by-layer deposition process, which is a feature of the fused deposition modeling (FDM) technique. Specifically, the tensile modulus and strength in the raster direction, the transverse direction, and the build direction are all evidently different [18,32]. To approach the reality, topology optimization for AM should deal with these directional variations. On the other hand, optimizing the filament deposition paths will provide an extra possibility to enhance the structural performance. However, the anisotropic material properties were often ignored by the topology optimization works for FDM printing [33]. Based on this background, a hybrid topology optimization method of design for multipatch FDM 3D printing is proposed in this paper, in which a printing layer is disintegrated into multiple patches and each patch has its unique direction of the zigzag-type filament deposition. In this way, the topology optimization of the targeted problem evolves from the traditional material/void interface design to a more sophisticated problem involving multiple levels of design freedom listed as follows: Design domain: The design domain defines the bounds of the material's distribution. By topology optimization, it will be divided into the material domain and the void.
Material domain: Originally, the material domain consists of homogeneous materials. Due to the different deposition directions, the material domain is further decomposed into several submaterial domains which are shown as the patches in Figure 1.

Sub-material domain (patch):
Each sub-material domain is characterized by a unique filament deposition direction which is to be optimized.  Figure 2 illustrates the three levels of design variables. A novel optimization problem is configured in this paper to simultaneously address all the design variables. Specially, the material domain is defined by a level set function and then optimized through level set topology optimization. The sub-material domains are not clearly distinguished from the beginning; instead, the candidate deposition directions are interpolated under the discrete material optimization (DMO) scheme with density variables. The sub-material domains will gradually emerge with the penalization of the density variables. Moreover, the deposition directions of the patches will be simultaneously optimized since the material constitutive model is deposition direction-dependent. More details will be provided in the next section.

Sub-material domain (patch):
Each sub-material domain is characterized by a unique filament deposition direction which is to be optimized. Figure 2 illustrates the three levels of design variables. A novel optimization problem is configured in this paper to simultaneously address all the design variables. Specially, the material domain is defined by a level set function and then optimized through level set topology optimization. The sub-material domains are not clearly distinguished from the beginning; instead, the candidate deposition directions are interpolated under the discrete material optimization (DMO) scheme with density variables. The sub-material domains will gradually emerge with the penalization of the density variables. Moreover, the deposition directions of the patches will be simultaneously optimized since the material constitutive model is deposition direction-dependent. More details will be provided in the next section.
For topology optimization with material anisotropy, the methods are conventionally categorized into two groups: topology optimization with discrete raster angles [34][35][36] and topology optimization with continuous raster angles [37][38][39][40]. The former employed the pointwise different raster angles which provides the largest design space; however, the discrete angle optimization result cannot be directly used for 3D printing, since it is non-trivial to derive equidistant continuous tool paths from the disorganized raster angles. The latter provides directly useful continuous deposition paths with equidistance; however, the design space is severely restricted by the contour-offset path pattern, i.e., all the deposition paths follow the shape of the structural boundary. Therefore, the motivation of this paper is to develop a new method that could achieve a balance between the manufacturability (with continuous deposition paths) and the allowable design space (with multiple patches). domain is defined by a level set function and then optimized through level set topology optimization. The sub-material domains are not clearly distinguished from the beginning; instead, the candidate deposition directions are interpolated under the discrete material optimization (DMO) scheme with density variables. The sub-material domains will gradually emerge with the penalization of the density variables. Moreover, the deposition directions of the patches will be simultaneously optimized since the material constitutive model is deposition direction-dependent. More details will be provided in the next section.  In the following section, the review of the AM induced directional material properties and their applications in engineering analysis and optimization are introduced. The multilevel optimization problem definition which involves the design domain modeling, the material domain modeling, and the sub-material domain modeling is presented in Section 3. Correspondingly, the solutions to the problem are revealed in Section 4. The case studies demonstrated in Section 5 prove the proposed method enhances the structural performance effectively. The conclusion of the contribution and the prospect of future work are made at the end.

Literature Review
Because of the anisotropic properties induced by the layered manufacturing process in AM, a great number of studies have been dedicated to studying the influence of anisotropic behaviors on design activities [41][42][43]. The relevant literature in this field is reviewed in this section.
Ahn et al. [32] investigated the anisotropic material properties of acrylonitrile butadiene styrene (ABS) parts made by FDM. Compared to the transverse direction, the material's tensile and compressible strengths are found to be stronger in the raster direction. Further, the classical lamination theory and Tsai-Wu failure criterion were applied to predict the failure load of FDM parts based on the experimental anisotropic material properties, which showed reasonable agreement with the experiments [44]. Lee et al. [45] conducted the compressive test of the FDM part built in different directions that observed evidently different compressive strengths. It has also been verified by the experiments conducted by Hill and Haghi [46] that the tensile modulus and tensile strength of the FDM polycarbonate are stronger in the raster direction.
Therefore, the build/raster directions should be taken care of to improve the part design due to the aforementioned anisotropic material properties. A novel cross-sectional structural analysis approach was introduced by Umetani and Schmidt [47], which emphasized the geometric relationship between the cross sections and the external load. Thus, the weakest cross section could be identified and the build direction is designed to be perpendicular to it. Ulu et al. [17] tested the anisotropic material properties through physical experiments. Then, surrogate modeling was performed based on finite element simulations of different build directions. Based on the constructed surrogate model, the optimal build direction was derived, which leads to a significantly improved factor of safety. It has also been proven that the derived structural performance can be significantly improved by simultaneously optimizing the build direction during topology optimization [19]. In a recent work, Mirzendehdel et al. [48] employed the Tsai-Wu failure criterion to constrain the anisotropic strength. The optimization result showed evidently more strength than the result design with the von Mises criterion. However, in this work, the build direction is not treated as a design variable.
Given raster direction optimization, topology optimization with discrete angle variables has attracted the major attention. Hoglund and Smith [34] performed concurrent structural topology and raster direction optimization, but the raster direction was treated as identical all over the design domain. Jiang et al. [35] later performed the similar concurrent optimization with pointwise different raster orientations. Yan et al. [36] developed a hybrid stress-strain method for the current bi-directional evolutionary structural optimization method (BESO) type topology optimization and the discrete raster angles, which addressed repeated global minimum problem. For the above studies, an obvious limitation is that the discrete angle optimization result cannot be directly used for 3D printing, since it is nontrivial to derive equidistant continuous tool paths. To overcome this limitation, topology optimization with continuous tool path planning was recently brought forward [35][36][37]. Most of the relevant studies focused on the contour-offset tool path which meant the filament deposition paths were created by offsetting the structural boundary. This type of tool path planning is made possible by employing the level set function for structural representation due to the signed distance feature [49]. However, the design space is severely reduced because the raster directions have to follow the shape of the boundary profile [43]. To avoid this issue, the multipatch tool path planning will be introduced in this paper. Specifically, the material domain will be divided into different patches and a unidirectional zigzag deposition path will be defined inside each patch. Then, the material domain representation and the patch (sub-material domain) representation are separated. As a result, a larger design space could be explored compared to design with contour-offset tool path, while the equidistance and continuity requirements can still be satisfied to guarantee the printability.

Design Domain Modeling
The level set function is applied to the model and optimize the material domain in this work. In [50], Osher and Sethian introduced the level set function as a natural way to represent a closed boundary. Specifically, if D is the initial design domain, Ω ∈ R n (n = 2 or 3) represent the material domain, and ∂Ω behaves as the material/void interface, then the level set function is defined as: The Heaviside and Dirac delta functions are commonly applied to model the domain and boundary, as shown in Equations (2) and (3) respectively: Then, the interior and boundary of the material domain can be expressed as: In practice, the approximate Heaviside and Dirac delta functions [51,52] are usually employed to facilitate the numerical domain and boundary integrations.
Level set also supports the multimaterial domain modeling. Multiple level set functions are available among which there are two generally applied approaches: the "color" level set [53,54] and the MMLS (Multimaterial level set) [55,56], which can be potentially extended to multimaterial multipatch FDM 3D printing.

Material Domain Modeling
As mentioned above, the material domain is decomposed into sub-domains in which the orientation of the zigzag-type filament deposition is distinct. In each sub-domain, by treating the material type and the corresponding raster direction together as an independent material type, multimaterial level set modeling makes it feasible to model the sub-domains. However, the authors intend to reserve the multimaterial level set modeling for the potential extension to multimaterial multipatch FDM. Therefore, another multimaterial modeling method DMO is applied in this work. DMO provides an effective way to conduct local interpolation between the candidate raster directions and it forces each local point to approach a specific raster direction at convergence. One benefit of combining the level set and DMO is that only n level set functions and m density variables are required to model n material types plus m raster directions.
Initially proposed by [57,58], DMO was utilized to interpolate between the multiple material selections applied to composite laminate design. Later, it has been actively applied to multilevel structural optimization problems [59]. The specific interpolation scheme is presented in Equation (6): in which ρ e 1 and ρ e 2 are local material volume ratios (i.e., relative densities) of material 1 and 2 which are distinguished by the different raster directions. D 1 and D 2 are the elasticity tensors corresponding to the two material types. p is the penalization term which is usually larger than 3 (i.e., p ≥ 3). The characteristic of the DMO scheme is that, the local material composition will converge to either (ρ e 1 = 1, ρ e 2 = 0) or (ρ e 1 = 0, ρ e 2 = 1). In summary, ρ e 1 and ρ e 2 are the density variables to distinguish the sub-material domains. Once reaching the convergence, (ρ e 1 = 1, ρ e 2 = 0) corresponds to D 1 , and (ρ e 1 = 0, ρ e 2 = 1) corresponds to D 2 . D 1 and D 2 are the material elasticity matrixes for the same orthotropic material model with different orientations.
In addition, it is trivial to extend the DMO to interpolate multiple candidate material types, as:

Sub-Material Domain Modeling
In each sub-material domain, the material properties will be affected by a raster direction variable θ. With this variable, the interpolation presented in Equation (6) is changed into: in which: and: It should be noted that, D 0 is a 2D elasticity tensor in case that the raster direction aligns with the x-axis and θ is counted in the counter-clockwise direction. Note that, the optimization is similar to multiscale topology optimization, wherein the second scale is treated as the raster directions. The constitutive model can be similarly derived from the homogenization theory.

The Overall Problem Definition
So far, the design freedom in each level has been specified. As presented in Equation (11), the overall problem is defined, which adopts the compliance minimization as the design objective. This problem definition is built under the assumption that the material domain is divided into two patches. (11) in which a(u, v, Φ, D) is the energy bilinear form and l(v, Φ) is the load linear form. u is the deformation vector, v is the test vector, and e(u) is the strain.
is the space of kinematically admissible displacement field. p is the body force and τ is the boundary traction force. V max is the upper bound of the material volume.
The solution to this problem will be presented in the next section.

Problem Solution
The optimization problem will be solved by the gradient based approach. The sensitivity analysis will be decomposed into three parts including (i) level set function, (ii) local densities, and (iii) raster directions.
For the first part (i), the sensitivity analysis of the compliance-minimization problem under the level set framework is quite mature. Therefore, the result is given directly as shown in Equation (12) without detailed derivation: where R is called shape gradient density, L is the Lagrangian, and λ is the Lagrange multiplier, which is updated with the augmented Lagrange multiplier method. Then, by following Equation (13), L is forced to change in the descent direction, as indicated by Equation (14): For part (ii) and (iii), the partial derivative on ρ e i and θ i are demonstrated in Equations (15) and (16), respectively: Therefore, the sensitivity results would be: Appl. Sci. 2020, 10, 943 (18) in which t is the time variable, and, Subsequently, the design update can be performed based on the sensitivity results. The level set function will be updated with the velocities of Equation (13) by solving the Hamilton-Jacobi equation. This method is well established [51,52]. The density variables and the direction variables can be updated with the steepest descent method based on the sensitivity information of Equations (17) and (18), since they are not involved in the constraints. Figure 3 illustrates the flow chart of the solution process. Note that, the three types of design variables are concurrently optimized in each iteration. It is also possible to employ the sequential optimization strategy [60], which however requires three times of finite element analysis of each iteration that lowers the efficiency.
Appl. Sci. 2020, 10, 943 7 of 18 to employ the sequential optimization strategy [60], which however requires three times of finite element analysis of each iteration that lowers the efficiency. It is worth noting that, the optimization problem is highly non-convex which easily traps the optimization result at a local minimum. The evolution of the three types of design variables can hardly be coordinated due to the different magnitudes of sensitivity results. A common situation is that, the density variables and the orientation variables evolve and form patches before the structure fully evolves to have clear-identified load paths. Consequently, misconverged local areas will appear. Figure 4 shows two examples that the raster directions inside the circles do not align with the longitudinal directions of the structural members, since the density and orientation variables inside the local areas converge before the structure reaching a clear load path. To alleviate this issue, an asynchronous starting strategy is proposed by blocking the sensitivity analysis and design update of the local densities for the initial N optimization iterations (indicated by the blue-colored dash lines in Figure 3). The purpose is to let the structural geometry evolve first to form clear load paths, and thereafter, to optimize the density variable to classify the patches. In this way, the local optimum issue as demonstrated in Figure 3 can be avoided.  It is worth noting that, the optimization problem is highly non-convex which easily traps the optimization result at a local minimum. The evolution of the three types of design variables can hardly be coordinated due to the different magnitudes of sensitivity results. A common situation is that, the density variables and the orientation variables evolve and form patches before the structure fully evolves to have clear-identified load paths. Consequently, misconverged local areas will appear. Figure 4 shows two examples that the raster directions inside the circles do not align with the longitudinal directions of the structural members, since the density and orientation variables inside the local areas converge before the structure reaching a clear load path. To alleviate this issue, an asynchronous starting strategy is proposed by blocking the sensitivity analysis and design update of the local densities for the initial N optimization iterations (indicated by the blue-colored dash lines Appl. Sci. 2020, 10, 943 8 of 18 in Figure 3). The purpose is to let the structural geometry evolve first to form clear load paths, and thereafter, to optimize the density variable to classify the patches. In this way, the local optimum issue as demonstrated in Figure 3 can be avoided.

Case Studies
In this section, a few numerical cases are studied to demonstrate the effectiveness of the proposed computational design method. In all the numerical examples, structured mesh with elements of unit size (1 by 1) is employed.

Cantilever problem
Firstly, the cantilever problem is investigated to minimize the structural compliance under the maximum material volume ratio of 0.5. The design domain has the size of 100mm by 50mm, which is discretized into 100 by 50 quadrilateral elements. The boundary conditions applied are shown in Figure 5(a). A point force of 1 kN is applied. The solid material with the Young's Modulus of 1.3GPa in the raster direction and 0.8 GPa in the transverse direction is assumed to be used. In addition, the Poisson's ratio is 0.4, and the shear modulus is 0.15 GPa. Figure 5(b) depicts the directions and the rotation angle of the raster direction, which is defined positively in the counter-clockwise direction.  As mentioned earlier, for a part made by AM, there can be three different levels of design freedom in total. To fully reveal their influences, two different optimization schemes will be tested with increased design complexities, which are: (1) Topology optimization with a fixed uniraster direction of 90°, 45° or 0°; (2) Topology optimization with two flexible raster directions starting from ±45°. The optimization results are demonstrated in Figure 6 and Table 1.

Case Studies
In this section, a few numerical cases are studied to demonstrate the effectiveness of the proposed computational design method. In all the numerical examples, structured mesh with elements of unit size (1 by 1) is employed.

Cantilever Problem
Firstly, the cantilever problem is investigated to minimize the structural compliance under the maximum material volume ratio of 0.5. The design domain has the size of 100mm by 50mm, which is discretized into 100 by 50 quadrilateral elements. The boundary conditions applied are shown in Figure 5a. A point force of 1 kN is applied. The solid material with the Young's Modulus of 1.3GPa in the raster direction and 0.8 GPa in the transverse direction is assumed to be used. In addition, the Poisson's ratio is 0.4, and the shear modulus is 0.15 GPa. Figure 5b depicts the directions and the rotation angle θ of the raster direction, which is defined positively in the counter-clockwise direction.

Case Studies
In this section, a few numerical cases are studied to demonstrate the effectiveness of the proposed computational design method. In all the numerical examples, structured mesh with elements of unit size (1 by 1) is employed.

Cantilever problem
Firstly, the cantilever problem is investigated to minimize the structural compliance under the maximum material volume ratio of 0.5. The design domain has the size of 100mm by 50mm, which is discretized into 100 by 50 quadrilateral elements. The boundary conditions applied are shown in Figure 5(a). A point force of 1 kN is applied. The solid material with the Young's Modulus of 1.3GPa in the raster direction and 0.8 GPa in the transverse direction is assumed to be used. In addition, the Poisson's ratio is 0.4, and the shear modulus is 0.15 GPa. Figure 5(b) depicts the directions and the rotation angle of the raster direction, which is defined positively in the counter-clockwise direction.  As mentioned earlier, for a part made by AM, there can be three different levels of design freedom in total. To fully reveal their influences, two different optimization schemes will be tested with increased design complexities, which are: (1) Topology optimization with a fixed uniraster direction of 90°, 45° or 0°; (2) Topology optimization with two flexible raster directions starting from ±45°. The optimization results are demonstrated in Figure 6 and Table 1. As mentioned earlier, for a part made by AM, there can be three different levels of design freedom in total. To fully reveal their influences, two different optimization schemes will be tested with increased design complexities, which are: (1) Topology optimization with a fixed uniraster direction of 90 • , 45 • or 0 • ; (2) Topology optimization with two flexible raster directions starting from ±45 • .
The optimization results are demonstrated in Figure 6 and Table 1.  Based on the optimization results, some interesting observations can be concluded as follows.
(1) The optimization result with the raster direction of 0° outperforms the ones with the raster Figure 6. Cantilever optimization results (the blue color shows the converged ρ e 1 = 1, the red color shows the converged ρ e 2 = 1, and the green color shows the misconverged or void areas). Based on the optimization results, some interesting observations can be concluded as follows.
(1) The optimization result with the raster direction of 0 • outperforms the ones with the raster directions of 45 • and 90 • . This result is reasonable because the principal stresses distribute along the horizontal direction of the beam in the presented cantilever bending problem.
(2) The optimization result with two designable raster directions yields the best structural performance. Specifically in Figure 6e, the area in blue color employs the raster direction of 15.44 • while the area in red color has the raster direction of −15.44 • .
(3) A clear-cut solid structural design is derived since level set method is employed. Then, within the solid area, ρ e 1 ≥ 0.95 or ρ e 2 ≥ 0.95 indicates clearly-formed patches, and 96.83% of the solid elements have the density variables properly converged, which validates the effectiveness of the DMO interpolation.
(4) The optimized shape in Figure 6c is non-symmetric about the horizontal axis even though symmetric boundary and loading conditions are applied. The reason for this non-symmetric phenomenon is that the material properties involved are not symmetric about the horizontal axis.
(5) Even though the raster directions are different, the results shown in Figure 6b-d share the identical topology only with some small variations in shape. Theoretically, topology changes occur by merging the predefined holes in the input topology structure shown in Figure 6a. Because the finite element discretization is 100 × 50, a dense interior hole distribution is not adopted considering the numerical stability. Therefore, the possible topology variations are restricted. On the other hand, the impacts of the different raster directions are reflected in the shape variations.

Short Cantilever Problem
Next, the optimization of short cantilever problem is conducted. The boundary conditions applied are presented in Figure 7a. A point force of 1 kN is applied. The problem configuration and the material properties are defined the same as the previous example, except that the maximum material volume ratio is 0.25 in this case. (4) The optimized shape in Figure 6(c) is non-symmetric about the horizontal axis even though symmetric boundary and loading conditions are applied. The reason for this non-symmetric phenomenon is that the material properties involved are not symmetric about the horizontal axis.
(5) Even though the raster directions are different, the results shown in Figure 6(b-d) share the identical topology only with some small variations in shape. Theoretically, topology changes occur by merging the predefined holes in the input topology structure shown in Figure 6(a). Because the finite element discretization is 100 × 50, a dense interior hole distribution is not adopted considering the numerical stability. Therefore, the possible topology variations are restricted. On the other hand, the impacts of the different raster directions are reflected in the shape variations.

Short cantilever problem
Next, the optimization of short cantilever problem is conducted. The boundary conditions applied are presented in Figure 7(a). A point force of 1 kN is applied. The problem configuration and the material properties are defined the same as the previous example, except that the maximum material volume ratio is 0.25 in this case.  The two different optimization schemes implemented in the previous example are investigated again. The corresponding optimization results are demonstrated in Figure 8 and Table 2. The two different optimization schemes implemented in the previous example are investigated again. The corresponding optimization results are demonstrated in Figure 8 and Table 2.   Through the analysis of the optimization results, both similar and exclusive conclusions can be drawn as compared to the previous example.
(1) The optimization results with the fixed raster directions of 0°, 45°, and 90° have nearly the same structural performance.
(2) The optimization result with two designable raster directions has a much better structural performance.
(3) A clear-cut solid structural design is derived since the level set method is employed. Then, within the solid area, ≥ 0.95 or ≥ 0.95 indicates clearly-formed patches, and 97.39% of the solid elements have the density variables properly converged.
(4) As shown in Figure 8(c), the optimized shape is still non-symmetric about the horizontal axis due to the non-symmetric material properties.

Michell structure
Then, the Michell structure problem is studied. The boundary conditions applied are shown in Figure 9. A point force of 1 kN is applied. The problem setup and the material properties are identically defined as the short cantilever case. Through the analysis of the optimization results, both similar and exclusive conclusions can be drawn as compared to the previous example.
(1) The optimization results with the fixed raster directions of 0 • , 45 • , and 90 • have nearly the same structural performance.
(2) The optimization result with two designable raster directions has a much better structural performance.
(3) A clear-cut solid structural design is derived since the level set method is employed. Then, within the solid area, ρ e 1 ≥ 0.95 or ρ e 2 ≥ 0.95 indicates clearly-formed patches, and 97.39% of the solid elements have the density variables properly converged.
(4) As shown in Figure 8c, the optimized shape is still non-symmetric about the horizontal axis due to the non-symmetric material properties.

Michell Structure
Then, the Michell structure problem is studied. The boundary conditions applied are shown in Figure 9. A point force of 1 kN is applied. The problem setup and the material properties are identically defined as the short cantilever case. Moreover, this example will highlight a new scheme of topology optimization with three designable raster directions (starting from 0° and ±45°). Scheme 2 defined in the previous two examples will also be studied for comparison. Correspondingly, the optimization results are presented in Figure 10. Moreover, this example will highlight a new scheme of topology optimization with three designable raster directions (starting from 0 • and ±45 • ). Scheme 2 defined in the previous two examples will also be studied for comparison. Correspondingly, the optimization results are presented in Figure 10. Moreover, this example will highlight a new scheme of topology optimization with three designable raster directions (starting from 0° and ±45°). Scheme 2 defined in the previous two examples will also be studied for comparison. Correspondingly, the optimization results are presented in Figure 10.  In comparison, the optimization with three designable raster directions depicted in Figure 10(b) has the structural compliance 7.38% smaller than the one with two designable raster directions. To be specific, the optimized raster directions of the design in Figure 10(a) are ±56.91°, while the optimized raster directions of the design in Figure 10(b) are ±57.18° and 0°.

Messerschmidt-Bölkow-Blohm (MBB) structure
At last, the MBB structure problem is studied. The boundary conditions applied are shown in Figure 11. A point force of 1 kN is applied. The problem setup and the material properties are In comparison, the optimization with three designable raster directions depicted in Figure 10b has the structural compliance 7.38% smaller than the one with two designable raster directions. To be specific, the optimized raster directions of the design in Figure 10a are ±56.91 • , while the optimized raster directions of the design in Figure 10b are ±57.18 • and 0 • .

Messerschmidt-Bölkow-Blohm (MBB) Structure
At last, the MBB structure problem is studied. The boundary conditions applied are shown in Figure 11. A point force of 1 kN is applied. The problem setup and the material properties are identically defined as the previous examples, except that the maximum material volume ratio is 0.4 in this case.
Similar to the Michell case, this example will explore the topology optimization with three designable raster directions (starting from 0 • and ±45 • ). Figure 2 defined in the first two examples will also be studied for comparison. Correspondingly, the optimization results are presented in Figure 12.  Similar to the Michell case, this example will explore the topology optimization with three designable raster directions (starting from 0° and ±45°). Figure 2 defined in the first two examples will also be studied for comparison. Correspondingly, the optimization results are presented in Figure 12. In comparison, the optimization with three designable raster directions depicted in Figure 12(b) has the structural compliance 5.92% smaller than the one with two designable raster directions. To be specific, the optimized raster directions of the design in Figure 12(a) are 39.99° for the blue colored area and −5.03° for the red colored area, while the optimized raster directions of the design in Figure  12(b) are 39.99° for the blue colored area, −48.12° for the red colored area, and 1.23° for the yellow identically defined as the previous examples, except that the maximum material volume ratio is 0.4 in this case. Similar to the Michell case, this example will explore the topology optimization with three designable raster directions (starting from 0° and ±45°). Figure 2 defined in the first two examples will also be studied for comparison. Correspondingly, the optimization results are presented in Figure 12. In comparison, the optimization with three designable raster directions depicted in Figure 12(b) has the structural compliance 5.92% smaller than the one with two designable raster directions. To be specific, the optimized raster directions of the design in Figure 12(a) are 39.99° for the blue colored area and −5.03° for the red colored area, while the optimized raster directions of the design in Figure  12(b) are 39.99° for the blue colored area, −48.12° for the red colored area, and 1.23° for the yellow Figure 12. MBB structure optimization results (the blue color shows the converged ρ e 1 = 1, the red color shows the converged ρ e 2 = 1, the yellow color shows the converged ρ e 3 = 1, and the green color shows the misconverged or void areas).
In comparison, the optimization with three designable raster directions depicted in Figure 12b has the structural compliance 5.92% smaller than the one with two designable raster directions. To be specific, the optimized raster directions of the design in Figure 12a are 39.99 • for the blue colored area and −5.03 • for the red colored area, while the optimized raster directions of the design in Figure 12b are 39.99 • for the blue colored area, −48.12 • for the red colored area, and 1.23 • for the yellow colored area. Moreover, the computing time of the case with two designable raster directions is 221s while the case with three raster directions consumes 585s. A major reason of the longer time of the latter case is that, it takes more iterations to make the density variables converge since three set of densities are involved. A desktop computer with an Intel Core I5-7400 CPU and 8GB RAM is used.
At the end, we would add a discussion about the initial-guess dependency issue. A finer mesh of 150 by 75 elements is used to perform the same topology optimization study with two designable raster directions. Initial setups of the level set function, density variables and raster directions all remain the same. Then, the optimization results are demonstrated in Figure 13. The optimized raster directions are 40.02 • for the blue colored area and −4.84 • for the red colored area. By comparing Figures 12a and 13, we can see that both the structural geometry and the raster directions are very close. The reason is that, the same initial guess of hole distribution is adopted for the two optimization schemes, while for the level set method, the design result is more sensitive to the initial artificial hole distribution instead of the mesh scale. Even so, the optimization with a finer mesh takes 165 iterations to converge which consumes 10 more iterations than the case with a coarse mesh; and the structural compliance of the finer mesh optimization is 1.51% smaller. The computing time of the finer mesh optimization is 675s.
Appl. Sci. 2020, 10, 943 15 of 18 colored area. Moreover, the computing time of the case with two designable raster directions is 221s while the case with three raster directions consumes 585s. A major reason of the longer time of the latter case is that, it takes more iterations to make the density variables converge since three set of densities are involved. A desktop computer with an Intel Core I5-7400 CPU and 8GB RAM is used. At the end, we would add a discussion about the initial-guess dependency issue. A finer mesh of 150 by 75 elements is used to perform the same topology optimization study with two designable raster directions. Initial setups of the level set function, density variables and raster directions all remain the same. Then, the optimization results are demonstrated in Figure 13. The optimized raster directions are 40.02° for the blue colored area and −4.84° for the red colored area. By comparing Figure 12(a) and Figure 13, we can see that both the structural geometry and the raster directions are very close. The reason is that, the same initial guess of hole distribution is adopted for the two optimization schemes, while for the level set method, the design result is more sensitive to the initial artificial hole distribution instead of the mesh scale. Even so, the optimization with a finer mesh takes 165 iterations to converge which consumes 10 more iterations than the case with a coarse mesh; and the structural compliance of the finer mesh optimization is 1.51% smaller. The computing time of the finer mesh optimization is 675s.

Conclusion
A hybrid topology optimization technique is presented in this paper in which the level set and DMO approaches are simultaneously applied to design and optimize the multipatch FDM parts. With the proposed technique, the concurrent material domain optimization, material domain decomposition with distinguished raster directions, and multiraster angle optimization have been realized. The effectiveness has been proven by a few case studies. In particular, aided by the proposed asynchronous starting strategy, the local misconvergence phenomenon of material domain decomposition has been successfully eliminated. It has also been observed that, the more patches simultaneously involved, the better the design performance will be.
In the future, the optimally designed parts are expected to be printed out so that their structural performance can be validated through experimental tests. In addition, the potential extension to multi-material, multipatch additive manufacturing will be explored. In the current study, a unidirectional zigzag deposition path is defined inside each patch for the sake of simplicity. In fact, the deposition path pattern of each patch can be extended to more complex patterns other than the zig-zag to achieve an even better design performance. Therefore, this aspect will be explored in our future work as well.

Conclusions
A hybrid topology optimization technique is presented in this paper in which the level set and DMO approaches are simultaneously applied to design and optimize the multipatch FDM parts. With the proposed technique, the concurrent material domain optimization, material domain decomposition with distinguished raster directions, and multiraster angle optimization have been realized. The effectiveness has been proven by a few case studies. In particular, aided by the proposed asynchronous starting strategy, the local misconvergence phenomenon of material domain decomposition has been successfully eliminated. It has also been observed that, the more patches simultaneously involved, the better the design performance will be.
In the future, the optimally designed parts are expected to be printed out so that their structural performance can be validated through experimental tests. In addition, the potential extension to multi-material, multipatch additive manufacturing will be explored. In the current study, a unidirectional zigzag deposition path is defined inside each patch for the sake of simplicity. In fact, the deposition path pattern of each patch can be extended to more complex patterns other than the zig-zag to achieve an even better design performance. Therefore, this aspect will be explored in our future work as well.