|Computer Modeling in Engineering & Sciences|
Topology Optimization of Self-Supporting Structures for Additive Manufacturing with Adaptive Explicit Continuous Constraint
College of Safety Science and Engineering, Civil Aviation University of China, Tianjin, 300300, China
*Corresponding Author: Jun Zou. Email: email@example.com
Received: 04 November 2021; Accepted: 29 December 2021
Abstract: The integration of topology optimization (TO) and additive manufacturing (AM) technologies can create significant synergy benefits, while the lack of AM-friendly TO algorithms is a serious bottleneck for the application of TO in AM. In this paper, a TO method is proposed to design self-supporting structures with an explicit continuous self-supporting constraint, which can be adaptively activated and tightened during the optimization procedure. The TO procedure is suitable for various critical overhang angles (COA), which is integrated with build direction assignment to reduce performance loss. Besides, a triangular directional self-supporting constraint sensitivity filter is devised to promote the downward evolution of structures and maintain stability. Two numerical examples are presented; all the test cases have successfully converged and the optimized solutions demonstrate good manufacturability. In the meanwhile, a fully self-supporting design can be obtained with a slight cost in performance through combination with build direction assignment.
Keywords: Topology optimization; additive manufacturing; self-supporting constraint; build direction assignment; gradual evolution
Additive manufacturing (AM) is a free-form manufacturing technique that builds parts in layers by depositing fine powder material, which is advantageous in building parts with complex shapes without specific tooling or fixturing . However, the AM is not an economical technology for conventional design, and the advantages can only be brought into play adequately through the elaborate design of high-performance structures. Topology optimization (TO) is an advanced structural design method that can be adopted to find the optimal material distribution within a given design domain under certain constraints. The state-of-the-art TO methods include homogenization method , solid isotropic material with penalization (SIMP) [3,4], evolutionary structural optimization (ESO) , level set method , moving morphable components (MMC)  and moving morphable void (MMV) . The integration of TO and AM technologies shows the potential to create significant synergy benefits .
In general, the results produced by conventional TO methods may not be built directly via AM. There are some manufacturing limitations that should be catered, such as minimum feature size, minimum slot distance, overhang limitation, material anisotropy, and microstructure control [1,9]. The lack of AM-friendly TO algorithms is a serious bottleneck for the application of TO in AM . Liu et al.  and Zhu et al.  summarized the state-of-art researches conducted on the integration of TO and AM for a variety of topics in recent years.
Overhang constraint is one of the most troublesome manufacturing limitations that are inherent within AM. When the overhang portions of a structure have inclination angles smaller than the critical overhang angle (COA), auxiliary supports are required to provide mechanical support and/or promote heat dissipation, otherwise the part will collapse during the AM process. Thomas  identified 45° as the COA for the parts manufactured with 316L stainless steel powder using SLM, while Wang et al.  found out that the overhang surfaces can always be well-printed if the overhang angles exceed 34°. Leary et al.  adopted 40° as the COA for those parts manufactured by fused deposition modelling (FDM). The auxiliary supports need to be removed after manufacturing, which will extend the build time, and increase material cost. The ideal approach is to integrate the COA restriction in the TO process, so that the solutions are self-supported and eliminates the need for additional supports.
Within the SIMP framework, Serphos  integrated the geometrical restriction in the optimization process through three strategies: multiple objective functions, global explicit constraint and density filter. Gaynor et al. [15,16] presented a support region projection to ensure a feature that is adequately supported, Then, Johnson et al.  made some improvements to allow for specifying COAs through an updated support region. Langelaar [18,19] introduced a self-supporting filter to rigorously exclude the geometries that violate the overhang angle criteria. Zhao et al.  devised an explicit quadratic continuous constraint to represent the self-supporting requirement. Qian  developed a Heaviside projection integral form with explicit geometric meanings based on the density gradient. Then, Wang et al.  introduced an approach that can be taken to optimize the build direction and topological layout simultaneously. van de Ven et al. [23,24] introduced an overhang filter with which the overhang regions can be detected by an anisotropic speed function. Garaigordobil et al.  proposed an overhang constraint based on an edge detection algorithm developed in the field of image analysis and processing. Zhang et al.  developed a method that can be applied to estimate the structural boundary normal by fitting local element density distribution with linear surfaces. Luo et al.  controlled the minimal overhang angle only in the enclosed voids to reduce the performance loss, and the inward-pointing normal is evaluated by the density gradient.
In addition to the SIMP frameworks, Mirzendehdel et al.  put forward a topological sensitivity approach for constraining support structure volume with in the level set framework. Liu et al.  proposed a multi-level set interpolation to address the support-free manufacturability constraint. By adopting polygon-featured holes as the basic design primitives, Zhang et al.  introduced the overhang constraints directly into the geometry description. Wang et al.  proposed a new form of overhang constraint in the level set framework, which is expressed as a single domain integral instead of point-wise constraints. Based on the MMC and MMV frameworks, Guo et al.  addressed the overhang limitation by constraining the angle of the curve boundaries. Zhao et al.  proposed a density filter within homogenization theory-based method to ensure that the porous structure satisfies the self-supporting constraint. Wang et al.  presented a B-spline based method to design self-supporting structures, in which an overhang angle constraint and a triangle constraint are presented.
This article presents an extension and improvement of the work presented in Zou et al. , in which an explicit self-supporting constraint model is constructed. The non-trivial improvements made in this article are as follows. Firstly, the improved method is applicable to various values of COA other than 45°. Secondly, the TO optimization process is combined with build direction assignment to reduce the performance loss. Thirdly, an adaptive explicit continuous constraint is devised to control the stability of structural evolution. Lastly, a triangular directional self-supporting constraint sensitivity filter is proposed to promote the downward evolution of structures and maintain stability. Furthermore, the volume preserving Heaviside projection filter  is adopted to suppress the gray elements and ensure the stability of optimization.
The remainder of this article is organized as follows. Section 2 introduces the formulation of the adaptive explicit self-supporting constraint in the SIMP framework. Section 3 introduces the problem formulation and TO procedure combined with build direction assignment. In Section 4, two numerical examples are presented to demonstrate the effectiveness of the proposed method. The conclusions are drawn in Section 5.
2 Self-Supporting Constraint Formulation in SIMP Framework
2.1 Brief Introduction of SIMP
The SIMP method is a widely used TO method, in which the elements of the discretized design space are assigned with continuous variables such as material density. The material properties such as Young's modulus are a function of the element density given by the interpolation scheme. In general, a penalization factor is introduced by raising the density to the power of this penalization factor to ensure black-and-white solutions. Besides, filtering methods should be applied to avoid checkerboards and mesh-dependency solutions. Compliance minimization is one of the most common TO problems applied for weight saving in various industries, which can be formulated as Eq. (1):
where c represents the compliance, F indicates the force vector and U denotes the nodal displacement vector. and refer to the design variable and physical density, respectively. K denotes the global stiffness matrix. V and V0 are the material volume and design domain volume, respectively, and Vf is the prescribed volume fraction.
A SIMP interpolation scheme proposed in Sigmund  is introduced to ensure black-and-white solutions, as formulated in Eq. (2):
where E0 is the Young's modulus after the interpolation. Emin is the stiffness of void material to void singularity of the stiffness matrix, which is generally considered as Emin = 10−9. p denotes the penalization factor, which is taken as 3 generally.
2.2 Adaptive Self-Supporting Constraint Formulation
The continuous characteristics of the design domain, such as element density gradient, can be relied on to identify the structural boundary and inclined angles, while this additional step is expected to make the optimization problem overly complex. In this article, the discrete element density is adopted to identify the unsupported elements. This article only considers the 2D rectangular region problems, in which all four sides of the design domain can be considered as the potential baseplate. They are denoted respectively as B, R, U, and L, as shown in Fig. 1. In order to adapt to various COAs, the design domain is discretized into rectangular elements, in which the element aspect ratio b/a is made adjustable according to the magnitude of COA and build direction. Take the bottom side as an example. Each solid element should be sufficiently supported by the three adjacent elements in the underlying layer.
The mathematical model of self-supporting criterion is defined as:
where is the density of element (i, j) and is the maximum density of the supporting elements. A smooth P-norm approximation is performed to represent the maximum function, and a linear penalization is performed to reduce the approximation error , in which the parameter pn controls the approximate accuracy and smoothness.
Eq. (3) can be used to judge the supporting status of each elements, and then an explicit constraint function is devised as Eq. (4) to quantify the global degree of self-supporting constraint violation, which will be added to the TO formulation.
where parameter is adopted to exclude the interference caused by density gradual transition and approximation errors, and it should be small enough to recognize the differences of density between each layer. Ae is the element area, which is used to make the constraint values evaluated in different build directions comparable. Parameter k is used to control the functional steepness. It is worth noting that when k is large enough, and .
For a clear black-and-white topology, it can be seen clearly that for self-supported design. While in practical numerical applications, the intermediate density elements are unavoidable in SIMP, and a density gradual transition will exist between the solid and void regions caused by the density filters. In addition, due to the approximation error, the value of will be small but not zero even for a self-supported design. To deal with this problem, the constraint is slacked by a positive number , which is restated as follows:
Different from the compliance, the self-supporting constraint is a local feature. In the early optimization iterations, there exist a large number of intermediate elements, and the structural boundary is blurry. In this circumstance, the constraint evaluation is inaccurate and meaningless. In order to make the optimization process more stable and prevent the local optimal solution with poor performance, a reasonable idea is to make the upper bound of constraint assigned and tightened adaptively during the optimization process. Therefore, the value of ɛ is large enough first to keep the self-supporting constraint inactivated until the structural compliance becomes stable. Then, the constraint is activated artificially and the following scheme is used for ɛ to decrease from to gradually:
where loop0 and represent the number of iteration and the corresponding value of when self-supporting constraint is activated, respectively. loop is the iteration number, and is a parameter that is larger than 1 and controls the decreasing rate of ɛ. is a properly selected positive number that is less than 1, which will be studied in the numerical examples.
3 Numerical Implementation
3.1 Topology Optimization Formulation Combined with Build Direction Assignment
The formulation of minimum compliance structural topology optimization subject to volume constraint and self-supporting constraint is demonstrated in Eq. (7), which is combined with build direction assignment to reduce performance loss. The structural compliance, volume constraint and self-supporting constraint are evaluated based on the physical density .
where variable d is the build direction. For simplicity, only the 4 build directions, denoted as B, R, U, and L in Fig. 1, are taken into consideration. Nevertheless, test results show that the performance loss can be reduced in a large extent.
To avoid the formation of checkerboard patterns, intermediate field of densities is obtained using linear filter  as follows:
where Hei is the weight factor, Ne is the set of elements i for which the center-to-center distance (e, i) to element e is smaller than the filter radius R1.
Since the linear filter averages the density distribution and is not good to obtain an ideal black-and-white design, the volume preserving Heaviside projection filter proposed by Xu et al.  is applied after the linear filter in Eq. (8) to obtain the physical density .
where β is a parameter controls the smoothness of the approximation. The Heaviside function has such a property: if 0 ≤ < η, is equal to zero; if η < ≤ 1, is equal to one; if equals to η, is η. A continuation scheme is used where the parameter β is gradually increased during optimization to avoid local minima. The value of η is determined by bi-section search method to ensure the volumes before and after filtering are the same. This Heaviside projection filter in Eq. (9) can suppress the gray elements at the structural boundary, which will make the structural boundary clear while ensuring the stability in optimization.
In addition, it should be noted that the self-supporting constraint is a function of element densities and build direction d. When the constraint values are evaluated in different build directions, as shown in Fig. 1, the design domain should be remeshed according to the element aspect ratio and build direction, with the element densities reassigned by interpolation in line with the current design.
3.2 Sensitivity Analysis and Triangular Directional Sensitivity Filter
The optimization problem is solved by the gradient based algorithm Method of Moving Asymptotes (MMA) proposed by Svanberg . As for the sensitivity analysis of objective function c and volume constraint V, reference can be made to Andreassen et al. . The sensitivities of self-supporting constraint are derived as follows. Considering that the linear density filter and Heaviside projection operation are applied, the chain rule is adopted to compute the derivatives of with respect to the design variables.
The sensitivity analysis ignores the knock-on effect of unsupported elements, for which the change in will only affect the value of , and . The first term on the right side of Eq. (10) can be generally expressed as Eq. (11), in which each term can be computed according to Eqs. (3) and (5) by derivative product rule.
As for the second and third terms on the right side of Eq. (10), reference can be made to Xu et al.  and Andreassen et al. , respectively. It can be easily found out that the sensitivity values are nonpositive, which means that the constraint value is reduced by increasing the element densities below the unsupported elements, and self-supporting is realized by gradual downward evolution of structural boundaries.
In order to promote the downward evolution of structural boundaries and avoid tiny supporting structures, a directional self-supporting sensitivity filter is proposed in this article. Used in the previous works, a semicircular directional sensitivity filter will inherently lead to some unsupported elements at the bottom of semicircle, and cause disturbance to some degree. In this article, a triangular directional self-supporting constraint sensitivity filter is devised as follows:
where Ne is the set of elements within the triangular as shown in Fig. 2a. Hei is a weight factor. The unit of R2 and distance (e, i) is the number of elements, so that the actual shape of the triangle varies according to the element aspect ratio. As the example illustrated in Fig. 2b, the solid element is unsupported. Through the action of the proposed directional sensitivity filter, the modified sensitivities are all negative in the upside down triangle area, which promotes the downward evolution and avoids tiny supporting structures.
3.3 Topology Optimization Procedure
A flowchart of the TO optimization procedure with self-supporting constraint that is combined with build direction assignment is demonstrated in Fig. 3. Firstly, the value of ɛ is large enough to keep the self-supporting constraint inactivated, which means that the design space is optimized for minimum compliance in the early stage of optimization. If the change rate of compliance less than 0.01, the constraint values are evaluated in different build orientations as described in Section 3.1, and assign the build direction to the orientation that owns the minimum self-supporting constraint value. Then, the design domain is remeshed if necessary, and the element densities are reassigned by means of interpolation in line with the current design, and the self-supporting constraint is activated as described in Eq. (6). If the change rate of both compliance c and is less than 0.001 in 6 continuous iterations, the optimization process is terminated. Otherwise, and are updated according to Eq. (6) until the condition of convergence is satisfied.
4 Numerical Examples
Two numerical examples are discussed in this chapter to demonstrate the effectiveness of the developed optimization process. The material is isotropic with Young's modulus E0 = 1 and Poisson's ratio ν = 0.3. The self-supporting constraint parameters are set as k = 150, = 0.03. The parameter pn is set to 60. The initial relative densities are all set to 0.5. The move limit is set to 0.1. The width a of elements is fixed to unit length, and the height b of elements is adjustable according to the magnitude of COA and build direction. For simplicity, the filter radius R1 and R2 are made equal in all optimizations, so that the dimensions of filtering area are the same in width direction.
4.1 MBB Beam
The problem statements and boundary conditions of the MBB beam problem are illustrated in Fig. 4. The design domain is a rectangular area with a width of 150 and height of 50. A unit load is applied on the middle point of the top of the design domain. The initial value of Heaviside projection parameter β is set to 1 and doubled every 50 iterations.
For comparison, the MBB problem is first optimized without the self-supporting constraint. The volume restriction is taken to be 40% of the domain, and density filter is applied with a radius of 3.5. The optimal free design is depicted in Fig. 5 with compliance cref = 229.20, which is treated as a reference for the results with self-supporting constraint. It is easy to find out that the dome and bottom of the structures are totally horizontal, which cannot be manufactured in vertical direction. As a result of Heaviside projection, the transition zone becomes very narrow and leads to a clear structural configuration.
When the self-supporting constraint is considering, the COA is first set to 45°, while parameter and are set to 1.3 and 0.1, respectively. The initial build direction is in orientation B. As the self-supporting constraint is inactivated until the build direction is assigned, the initial build direction has little effect on the optimized result. Fig. 6 gives the design evolution history with self-supporting constraint, in which the green blocks represent the baseplate according to the build direction. Fig. 7 illustrates the convergence curves of compliance and self-supporting constraint value. Fig. 6g shows the optimum design after 235 iterations with = 2.15, and the final build direction is in orientation R. The minimum compliance c is 238.85, which is only 4.2% higher than the optimum result of free design.
To illustrate the evolution of structural boundaries, the main unsupported structures are indicated by red ellipses, and the moving directions are indicated by red arrows, as shown in Fig. 6. As can be seen from these results, the boundaries of the unsupported structures move gradually to the baseplate with the increase of iteration steps.
To illustrate the final result more clearly, the COA is marked by red dash lines. It can be noted that the self-supporting constraint becomes active in the final design, which indicates the validity of the proposed method.
The iteration curve of compliance is quite stable except for the fluctuations caused by the periodic increase of the Heaviside parameter β. Nevertheless, the fluctuations are very small by adopting the volume preserving Heaviside filter. However, the iteration curve of self-supporting constraint fluctuates in a large scale, which decreases significantly at iteration 40. The main reason for this is that the build direction has been assigned and self-supporting constraint is reevaluated. Besides, as the self-supporting constraint is imposed on the gradual density transition region at the structural boundary, which leads to a certain degree of fluctuations with the evolution of structural configuration. Nevertheless, the fluctuations are insignificant in the later period of optimization, which is mainly attributed to the triangular directional sensitivity filter. As the structural configuration eventually becomes clear, the constraint values become stable as well.
To investigate the effect of parameter δ on the solutions, the MBB beam is optimized with δ of 0.3, 0.5, 0.7, and 0.9, respectively. The values of all other parameters are the same as before. All optimizations are converged, and the final build directions are all in orientation R. For comparison, the solutions are rotated 90° clockwise, as shown in Fig. 8, where the elements in red indicate the structural boundaries that cannot be successfully printed. As suggested by these results, the constraint value increased almost linearly with the increase of δ. In fact, as the elements in red cannot be successfully printed, the adjacent solid elements above the red elements cannot be manufactured successfully either. As described in Section 2.2, the parameter δ controls the final upper bound of self-supporting constraint, and 0.1 is used in the following cases.
In order to fully demonstrate the performance and robustness of the proposed algorithm, the optimizations with varying filter radius and COA are studied. In general, the optimization with large COA is hard to converge as the manufacturing constraint becomes very stringent, while high overhang angles are usually easy to print. For simplicity but without the loss of generality, a parameter study for 30°–60° is performed. The MBB beam is solved with filter radius R of 3.5 and 4.5 and COA of 30°, 45°, and 60°. The performance of the proposed method in maintaining the structural performance is measured by compliance ratio, c/cref, which is the structural compliance computed with or without self-supporting constraint.
All optimizations converged and realized self-supported successfully, with the optimized solutions and performances shown in Fig. 9, in which the baseplate is also presented by green blocks. It can be seen that the final build direction varies with different COAs, which is natural as the COA has direct influence on the constraint value . Besides, it is worth noting that the build directions are different for filter radius of 3.5 and 4.5 when COA = 30°. In addition, the results suggest that a smaller filter radius will result in finer structures, which is consistent with what has been seen in any free-form topology optimization.
As can be seen in Fig. 9, the compliance ratio c/cref is much smaller than the results in , which validates the effectiveness of the proposed TO procedure combined with build direction assignment. As expected, the result obtained when COA = 30° and filter radius is 3.5 is almost the same as the freely optimized result shown in Fig. 5. With the increase of COA, the constraint becomes stricter and makes the solutions deviate significantly from the free design results, which leads to greater performance loss. In addition, it can be noticed that the elements size with COA = 30° is larger than that of 60°. The reason for this is that the values of element height b are different. As the width a of elements is fixed to unit length, the elements height b is for COA = 30° when the build direction changes to horizontal direction, while it is only for COA = 60°.
4.2 Cantilever Beam
The second case is a cantilever beam problem with two subcases, as shown in Fig. 10. The design domain is a 150 × 50 rectangular area, while in the general domain at right there is a prescribed circular hole at the center of design domain. Both models are fixed on the left edge with a unit external force exerted on the middle point of the right edge. The initial value of Heaviside projection parameter β is set to 0.25 and doubled every 50 iterations.
The two cases are both first optimized without the self-supporting constraint. The allowable volume fraction is set to 30%, and filter radius is set to 3.5. The results are depicted in Fig. 11, and the compliance is 275.74 and 307.89, respectively. The prescribed hole will provide a challenge for evolving self-supported structures, which is conducive to testing the weaknesses of the proposed method.
When the self-supporting constraint is considering, the COA is first set to 60°, while parameters and are set to 1.1 and 0.1, respectively. As mentioned before, a smaller value 0.25 is adopted as the initial Heaviside parameter β. The is because a smaller β will lead to a wider gradual density transition region, which is conducive to the gradual evolution of structural boundaries. The optimization of rectangular domain case converges at iteration 270 with = 4.73, and the final build direction is in orientation L. The optimal compliance c is 310.4 with a compliance ratio of 111.0%. Fig. 12 illustrates the design evolution history, and Fig. 13 presents the convergence curves of compliance and self-supporting constraint value.
Similarly, the evolution of structural boundaries is indicated by red ellipses and red arrows, while the COA is marked by red dash lines. As can be seen from these results, the unsupported structures move toward to the baseplate gradually until the structure becomes fully self-supported. Compared with the reference design, the number of staggered structures is smaller and the structural layout is clean and simple.
It can be noted in Fig. 13 that the fluctuating range of compliance curve is much higher compared with the MBB example. The reason is that the structure configuration changes a lot as compared to the free design results. Besides, there are many intermediate density elements created during the evolution of structural boundaries due to the move limit. The build direction is assigned at iteration 47, which successively results in a significant reduction in self-supporting constraint value. It can be found out that the structural layout at iteration 40 is almost the same as the optimum design in Fig 11a except that there are more intermediate density elements.
Like the first example, the rectangular domain case is also optimized with δ of 0.3, 0.5, 0.7, and 0.9, respectively. All optimizations converged, and the print direction are all in orientation R. Similarly, all solutions are rotated 90° clockwise, as shown in Fig. 14. According to this figure, with the increase of δ, the number of unsupported elements increases and the constraint value increases correspondingly. In the following cases, the parameter δ is set to 0.1.
In order to fully demonstrate the performance and robustness of the proposed algorithm, the rectangular domain case is also solved with a filter radius of 2.5 and 3.5 and a COA of 30°, 45°, and 60°. All optimizations realize self-supported successfully, with the results shown in Fig. 15. It can be seen that all build directions are in in orientation L. With the increase of COA, the structural configuration changes more significantly compared with the free design result, and lead to a more significant performance loss. As expected, the result of COA = 30° is most similar to the freely optimized result. The influence of filter radius on structure configuration is less obvious than the first example. Nevertheless, a larger number of intermediate density elements can be noticed at structure boundaries for the filter radius of 3.5.
In addition, the cantilever beam with non-designable hole is solved with a COA of 30°, 45°, and 60° and an allowable volume fraction of 30% and 40%. The values of all other parameters are the same as before. The solutions are shown in Fig. 16. It can be seen from the figure that all solutions are self-supported, and the final build direction is in L orientation. With the increase of volume fraction, the structural compliance c decreases in a large extent as more elements can be used for load transfer. Besides, like in the previous cases, the structural configurations change more significantly compared with the free design results for larger COA, which leads to greater performance loss. The results demonstrate that the proposed method is effective in dealing with the structures containing non-designable domains.
From the results as mentioned above, it can be found out that there are barely tiny or distorted structures in the solutions with self-supporting constraint. This is mainly attributed to the adaptive explicit continuous constraint, which realizes self-supported through the gradual evolution of supporting structures, and self-supported is only required at the final convergence point. The structural layouts are natural and clean, which has a significant meaning in improving manufacturability and alleviating stress concentration.
In this article, a TO procedure of self-supporting structures that suitable for various COA is proposed, which is integrated with build direction assignment to reduce performance loss caused by self-supporting constraint. Besides, a new adaptive explicit continuous self-supporting constraint function is constructed to quantify the degree of self-supporting constraint violation, which can be adaptively activated and tightened during the optimization process. In order to apply it to different COAs and build directions, the design domain is discretized with rectangular elements in a way that the element aspect ratio can be automatically altered. Besides, the volume preserving Heaviside projection filter is adopted to improve stability, and a triangular directional self-supporting constraint sensitivity filter is devised to promote the evolution of supporting structures and maintain stability.
Two numerical examples are presented to demonstrate the functionality and performance of the developed optimization procedure. All of the test cases have successfully converged and achieved self-supported, and the optimized solutions show good manufacturability. With the increase of COA, results show that the structural configuration changes more significantly compared with the result of free design, which leads to more severe performance loss. Nevertheless, by combining the TO optimization process with build direction assignment, the performance loss is reduced in a large extent.
In the future works, we will focus on extending the proposed method to 3D problems, and one of the main challenges is the computational efficiency. Besides, the definition of self-supporting criterion can be improved for unstructured meshes, so that the topology optimization can be performed in arbitrary build direction, and address complicated design cases.
Acknowledgement: The authors thank Krister Svanberg for the use of MMA optimizer.
Funding Statement: This work is supported by the National Key Research and Development Program of China (2018YFB1106303), and Scientific Research Foundation of CAUC (2017QD10S).
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|