An Accelerating Convergence Rate Method for Moving Morphable Components

In the structural topology optimization approaches, the Moving Morphable Components (MMC) is a new method to obtain the optimized structural topologies by optimizing shapes, sizes, and locations of components. However, the size of the mesh has a strong influence on the rate of which the component builds the initial topological configuration by moving. -e influence may slow down the convergence rate. In this paper, a hierarchical mesh subdivision solution method that can accelerate the convergence rate for the MMC is developed. First, the coarse mesh is used as the starting point for the optimization problem, and the construction process of the initial topology structure is increased speed by accelerating the movement of components. Second, the optimized solution obtained by the coarse mesh is equivalently mapped to the same problem with a finer mesh and used to construct a good starting point for the next optimization. Finally, two-dimensional (2D) MBB beam example and three-dimensional (3D) short cantilever beam example are provided so as to validate that with the use of the proposed approach, demonstrating that this method can improve the convergence rate and the stability of the MMC method.


Introduction
In general, structural topology optimization aims at finding the optimal distribution of the prescribed material volume within a given design area and boundary conditions to obtain the structural layout that meets the desired characteristics. Compared with size and shape optimization, structural topology optimization has more flexibility and is more challenging in dealing with the optimization problems. Since the pioneering work of Bendsøe and Kikuchi [1], structural topology optimization has attracted the attention of more and more researchers. After decades of development and improvement, many approaches have been advanced for structural topological optimization and some of them have successfully solved different problems in various physics disciplines such as acoustic [2], electromagnetic [3], and optics [4]. Nowadays, structural topology optimization has become a vital part of the engineering design process (for a state-of-the art review of structural topology optimization, see [5][6][7][8][9]).
Although existing topology optimization methods have the potential to provide novel concept designs with little prior knowledge, it is worth noting that from the geometry expression point of view, most of them are developed within the pixel (i.e., the presence or absence of density within finite elements) or node point-based solution framework. Among these approaches, the solid isotropic material with the penalization (SIMP) method [10][11][12] and level set method [13][14][15] are the most well-established methods. In the SIMP method, the design domain is first discretized into a number of finite elements (pixels) and the density is assigned to each element as a design variable. en, the optimal topology can be expressed by the presence or absence of density in pixels that are calculated using mathematical programming or optimality criteria algorithms. e optimized structures obtained by the SIMP method are prone to boundaries with gray-scale diffusion or jagged shapes. Compared with the SIMP method, in the level set method, the values of the level set function at the node points of the finite element mesh are often taken as the design variables, and the structural boundaries can be specified by drawing the contour of the level set function in one higher dimensional space. erefore, the level set often obtains the optimization structure with smooth boundaries [16]. However, whether the SIMP method or level set method, they all solve topology optimization problems with an implicit way (i.e., no accurate geometry information is embedded in them), which will make it difficult to accurately control the structural feature sizes.
In the MMC method, the precise geometric information can be optimized according to the structural responses during the optimization process. However, the optimization process has some convergence problems [35]. Recently, Norato et al. [36] applied the same idea of the MMC method to discrete elements by using projection techniques. When analyzing the dependence of the initial layout, they found that the convergence rate for different initial layouts is different. Zhang et al. [37] studied the influence of the initial component quantity and layout on the final structural complexity. ey found that the unreasonable positioning of the components may cause the oscillatory behavior of the convergence process, which will affect the rate of the convergence. Deng and Chen [38] used the connected morphable components (CMC) method in the initial layout to enhance the connectivity of the topology structure, and the stable and efficient convergence history can be ensured.
Although the above scholars have found that the MMC method has some convergence problems, most of them rarely do further research to improve the rate of the convergence. erefore, in order to improve the convergence rate of the MMC method, we propose a hierarchical mesh subdivision method.
is approach starts from a coarse finite element mesh, and the obtained optimal solution by the coarse mesh is used to construct a good starting point for the next finer mesh. is method can not only reduce the finite element calculation consumption but also accelerate the determination of the topology structure. Compared with the original MMC method, the convergence rate of the original MMC method is improved. e effectiveness of the proposed method is verified by 2D and 3D numerical examples.

MMC Method
In the MMC method, the structure can be constructed by the corresponding topology description function (TDF, φ(x)) of components in the following form: where R denotes the whole design domain and Ω w ⊂ R represents the domain occupied by all components containing solid material, i � 1, . . . , N indicates the number of the components, the TDF of the ith component can be express as where Ω w i denotes the design domain comprised by the ith component and Ω w � ∪ N i�1 Ω w i , and the geometric representation of the structure is shown in Figure 2. It is worth noting that the not differentiability of equation (2) is neglected in the original MMC method [18].
According to [17,18], the geometry of a component is represented by constructing the TDF. e half-length, the half-width, the coordinate of the center, and the tilt angle of components are considered as the design variables to control the components' movement and deformation. e TDF of the ith component can be denoted as where in which x i0 , y i0 , l i , sin θ, t i1 , t i2 , and t i3 are defined as in Figure 3. For the description of a 3D structural component, see [21].
In the study of this paper, the optimization problem is designed considering minimization of the structural compliance under a constraint on the available solid material. Based on the MMC method, the optimization problem can be written as where J denotes an objective function to be minimized, K, U, and F denote the global stiffness matrix, displacement vector, and load vector, respectively, D denotes the global design variable vector, D i , i � 1, 2, . . . , N represents the vector of design variables associated with the ith component, where N is the total number of components, and V max denotes the upper limit value of V(D).

A Hierarchical Mesh Subdivision Solution Method
In the MMC method, the lengths, the widths, the central coordinates, and the tilt angles of components are treated as design variables; thus, the structural topology change can be achieved by the movement, deformation, overlap, and merging of components within the design domain. Compared with the traditional topology optimization methods, the computational cost of variables is significantly reduced. However, the movement of the component in the MMC is achieved by the transformation of the central coordinate. As shown in Figure 4, the center coordinates of components are located on the mesh nodes with distance units, which are formed by mapping the original mesh nodes in a certain way. erefore, the rate of which the topology structure is constructed by component movement is related to the fineness of the mesh.
As shown in Figure 5, a design area with the same boundary conditions are discretized into different meshes, and the initial layout of the components is (2, 1, nelx, nely, 0.5, 0.5, [0.38, 0.08, 0.08, 0.08, 0.7], 0.4). e components with the black boundary are the initial layout and the components with red and blue boundaries are formed after two iterations of the adjacent, respectively. We take any two components I and J in the design domain, and (x, y)denotes the center coordinate of one component. en, the center coordinate position of the adjacent iteration step component can be expressed as where j represents the iteration step. e moving distances of the components I and J can be calculated by We take j � 1, 2; the relevant data is shown in Table 1.
It can be clearly seen from Table 1 that, in the same design domain, the finer the mesh is, the more nodes are included, but the moving speed of the components also may be slowed down under the finer mesh. Although the small distance movement of the components can increase the accuracy of the connection between components, in the initial stage this will not only increase the computational cost, but it may also slow down the rate of the convergence due to the difficulty in quickly determining the location of some components. If the coarse mesh is used to accelerate the movement of the components in the initial stage of   optimization. en, the obtained optimal solution is used to construct a starting point for the next finer mesh, which will speed up the convergence rate. erefore, a natural idea is to employ the multigrid approaches to the MMC framework. Several multigrid methods have been used in density-based topology optimization. e uniform element size grid strategy [39,40] applies the solution of the coarse mesh optimization to subsequent suboptimizations with the finer uniform element size grid. Adaptive Mesh Refinement (AMR) strategy [41][42][43] aims to have a coarser mesh in the void region, and a finer mesh for regions that are solid or near the boundary, which is more effective in saving computational consumption. As shown in Figure 6(a), in order to combine the AMR strategy with the MMC framework, finer grid cells are used at the boundary and inner areas of the components and coarser grid cells are used in the rest. In this case, components that need to move a small distance can meet their precise positioning within a finer grid. But when the distance in which some components need to move is between the finer grid node and the coarse grid node, because there are no suitable nodes inside the coarse grid, it may move to the closer or farther nodes (see Figure 6(b)).
is may affect the precise positioning of some components and slow down the rate of convergence. In addition, in the MMC method, only components whose thickness is greater than a mesh element size can be displayed. Just as shown in Figure 5(a), the boundary of the component has been deformed under the coarser mesh. In the AMF method, some components located in a finer grid are likely to be smaller than the coarser grid elements existing in the design domain. If they are moved from the finer grid area to the coarser grid area, it will cause a temporary missing display of some components in the structure topology, as shown in Figure 6(b). In comparison, the uniform distribution of nodes in the uniform element size grid method is more advantageous in the positioning and display of components. erefore, this paper applies a uniform size hierarchical meshing strategy for the MMC framework.   Based on the above analysis, a method of combining a uniform size grid strategy and the MMC framework is proposed. e main idea of this method is to divide the optimization process of the MMC method into several substages; this means that the sequence of problems, P q , that we need to solve can be described as where n q � 4n q − 1, n represents a positive integer, q denotes the number of mesh subdivisions, i.e., q � 0 indicates that the number of mesh is the initial value nelx × nely. Each subproblem P q is solved by the method proposed in this paper, and the obtained optimization solution of P q is used to construct a starting point for the next subproblem P q+1 . As shown in Figure 7, firstly, the design domain is first divided into nelx × nely mesh, and we assign the initial values of the design variables ini_val � [x 0 , y 0 , l, t 1 , t 2 , t 3 , sin θ] to the variable set of all components: Secondly, when the mesh subdivision condition is satisfied, the coarse mesh nelx × nely will be discretized into a finer mesh 2 × nelx × 2 × nely, and the topological solution optimized under the coarse mesh will be equivalently mapped to the same problem with a finer mesh to optimize. Finally, if the convergence condition is met, the final topology is the output; otherwise, the mesh subdivision and the optimization are continued.
In the hierarchical mesh subdivision method, the number of discrete mesh can be expressed as where nelx × nely denotes the number of meshes before subdivision, 2 × nelx × 2 × nely represents the number of meshes after subdivision, and Iters means the number of iterations. Frc denotes the mesh subdivision condition, which can be written as where Max iter denotes the maximum number of iterations, n j denotes a positive integer, and the reader can choose different values according to your needs, and j ∈ N represents the number of the mesh refinements. In the optimization process, when the mesh refinement condition is satisfied, each square element is divided into four new square elements. e number of meshes will be divided into 2 j × nelx × 2 j × nely mesh. As shown in Figure 8, the center coordinates of a component A ik in the structure obtained by the coarse mesh optimization are on the node (i, k). When the mesh is subdivided, each unit is subdivided into 4 units of 4 nodes, but the position of the nodes does not change after the mesh is subdivided. erefore, after the coarse mesh is subdivided, the equivalent transformation of the topology structure can be realized by mapping the components variables set Τ of the coarse mesh output to the node corresponding to the fine mesh. In the 3D case, the number of subdivided coarse meshes nelx × nely × nelz will be subdivided into a number of meshes 2 j × nelx × 2 j × nely × 2 j × nelz, and the same component variables equivalent mapping method can be used to achieve hierarchical mesh subdivision optimization. e biggest difference of the strategy is that the mesh can be subdivided as needed during the optimization process.

Finite Element Analysis.
In the paper, the Eulerian mesh is utilized and TDF is used to calculate the value of each node. e ersatz material model will be used to calculate the values of TDF at four nodes of an element, and Young's modulus of every element can be obtained by (11), as shown in [18]: where q is an integer number, we take 2, E 0 denotes Young's modulus for the material of each element, H � H(x) is a regularized Heaviside function, and φ e j , j � 1, . . . , 4, are the values of the TDF function of the whole structure at four nodes of an element: where ε is the regularization parameter and a � 0.001 to prevent the occurrence of the singular of the global stiffness matrix, and the stiffness matrix of the ith element can be written as where k s is the element stiffness matrix when E i � 1.

Sensitivity Analysis.
Once the stiffness matrix of each element is been obtained, the global stiffness matrix K can be assembled. e objective function C and the displacement vector U will be constructed as Assume that μ is any one of the variables, and its sensitivity of the structural compliance can be written as Since the force is independent of the design variables, we can get (17) by (11) and (13): By equations (15), (16), and (17), the derivative of compliance could be written as For the derivative of volume constraint, we also have

Numerical Examples
In this section, two numerical examples are investigated to demonstrate the effectiveness of the proposed method. Unless otherwise stated, all involved quantities are chosen as dimensionless. Uniform four node bilinear square elements are used to discretize the design domain. Young's modulus of the solid material is E � 1 and Poisson's ratio is υ � 0.3. e method of moving asymptotes is used as the numerical optimizer; the hyperparameters of MMA are epsimin � 10 − 10 , raa0 � 0.01, albefa � 0.4, asyinit � 0.1, asyincr � 0.8, and asydecr � 0.6. e computer code analysis is performed on the LENOVO-TY5060 workstation with an Intel(R) Xeon(R) i5-6400 2.70 GHz CPU, 4 GB RAM of memory, Windows 10 and run with MATLAB R2017a.

e 2D MBB Beam Example.
In order to demonstrate the effectiveness of the hierarchical mesh subdivision method in improving the convergence rate, the 2D MBB beam example is considered. e design domain, boundary conditions, and external load are shown in Figure 9. A unit vertical load is imposed at the middle point of the upper surface. e upper bound of available solid material volume constraint is V � 0.4. e 2D MBB beam problem is optimized by the MMC method, the MMC method with the hierarchical mesh refinement (MMC-HMR), and the SIMP method, respectively. e initial parameters of each component are all set as 0.38 0.04 0.06 0.04 0.7 , and the initial components layout is shown in Figure 10. Because of the natural symmetric property of the problem, only half of the structure is used to optimize. When using the MMC method and the SIMP method, the design domain is discretized by 240 × 80 mesh. In the MMC-HMR method, the design domain is first discretized into 60 × 20 coarse mesh, and the coarse mesh is refined 1 and 2 times, respectively. e sizes of the refined meshes that are actually involved in analysis are 120 × 40 and 240 × 80 mesh. In (10), j � 0, 1, 2, n 0 � 10, n 1 � 10/3, and n 2 � 5/3. In the SIMP method, the penalty parameter p � 3 and the filtering radius r � 6.
Some intermediate steps and the history of the convergence process of the three methods are presented in Tables 2 and 3. As shown in Tables 2 and 3, the SIMP method gradually optimizes the structural topology by penalizing the density of the elements. It obtains the structural topology with fewer iterations, and the entire optimization process is stable. However, its disadvantages are that the final structure does not contain accurate geometric information, and grayscale diffusion often exists at the boundary of the structure, which makes subsequent optimization and manufacturing of the structure difficult. Compared with the SIMP method, the MMC method and the MMC-HMR approach have clearer and smoother structural boundaries. From Table 2, we can clearly see that the structural topology of the MMC method has undergone major changes due to the unreasonable positioning of some components. is may not only slow down the convergence of the optimization process, but also affect the final structure. e history of the convergence process shows that the value of the objective function by using the MMC method has undergone oscillation erratically from 250 iterative steps to the maximum number of iterative steps. In contrast, the initial configuration of the structure can be quickly established with the 60 × 20 coarse FE mesh by the proposed MMC-HMR method, and the solution is used as a starting guess for the same problem with a refined mesh. e topology structure rarely undergoes significant mutations after this solution, which also enhances the stability of the optimization processes at next stages. e final optimal results are shown in Table 4. It is worth noting that, the structure obtained by the MMC method is slightly different from the result of the SIMP method, which may be due to the optimization process which does not converge to a stable state. e result obtained by the paper method is very similar to the result of SIMP. Although SIMP has the fewest iterations, it also costs the most time because each finite element contains a design variable. In the SIMP method, 102 iterations need 17645.88 s. e MMC method does not reach convergence within the set maximum number of iterations, and the total time cost is 2045.26 s. e MMC-HMR method reached convergence at the 639th iterative step. e total time cost of the paper method is 709.09 s, which saved about 65.3% compared with the time cost of the MMC method and saved about 95.98% compared with the time cost of the SIMP method.

A 3D Short Cantilever Beam Example.
e classical 3D short cantilever beam example is examined. e boundary conditions and external load are shown in Figure 11. Totally eighteen components are placed in the design domain. is problem is optimized by the MMC method, the MMC-HMR method, and the SIMP approach, respectively. e 3D MATLAB code of the SIMP method in [44] is used in this example. In the SIMP method, the penalty parameter p � 3 and the filtering radius r � 2.4. e initial components layout is shown in Figure 12.
e design domain is  e final optimal structural topology obtained by the different methods is shown in Figure 13. It can be clearly seen that some obvious fine components that are not connected and not removed in the results of the MMC method (see Figure 13(a)). In comparison, the structure obtained by the MMC-HMR method is very similar to that  Mathematical Problems in Engineering   of the SIMP method, and it is worth noting that the boundary of the optimal structure by the MMC-HMR is very smooth. Figures 14-16 show some intermediate results obtained during the optimization process of the three methods, respectively. e figures clearly show that the topology structure has undergone some drastic changes during the optimization process of the MMC method. In comparison, the topology structure obtained by the MMC-HMR method through the coarse mesh has hardly changed in the subsequent optimization process. e intermediate structure topology obtained by the SIMP method also has not changed much, but the disadvantages are similar to the previous 2D optimization result. e structure boundary is nonsmooth and contains more gray elements. Table 5 lists the value of the final objective function, the total number of iteration steps, the time cost, and the convergence history. To demonstrate the convergence process, the maximum iteration step of the MMC is set to  Figure 11: A 3D short cantilever beam example.  Figure 12: e initial design of the 3D short cantilever beam example. 4000. It can be clearly seen from Table 5 that the value of the objective function obtained by the SIMP method is much larger than that of the other two methods. Compared with MMC-HMR, because of the existence of more small structures in the structure, the MMC method achieved a slightly lower the value of the objective function. From the convergence history, we can see that the value of the objective functional obtained by MMC has oscillation behavior in the optimization process, and the optimization does not converge within the set maximum iteration step, which may be due to the uncertainty in the positioning of some components. In comparison, the convergence history of the SIMP method and the MMC-HMR method is more stable and requires fewer iterations. e SIMP method used only 100 iteration steps in this example. But its design variable is 10800, which is much larger than 162 in the MMC-HMR method, so the convergence rate is also slower than this paper method.

Conclusions
In the present paper, a hierarchical mesh subdivision method based on the MMC framework is proposed. Compared with the original MMC method, the proposed method divides the optimization process into several substages. e coarse mesh is used to accelerate the initial topological configuration, and the optimization result of the coarse mesh is utilized as a starting point for the same problem with a finer mesh. is not only accelerates the MMC optimization process but also enhances the stability. From the numerical examples of 2D MBB beam and 3D cantilever beam, compared with the original MMC method, the MMC-HMR method has a substantial improvement in convergence speed and stability. Although the proposed method is adopted on the MMC method, in principle other solution frameworks based on an explicit geometric description can also use the method to improve the convergence rate. In this paper, the conditions of mesh subdivision are set according to the number of iterations. A more efficient and accurate automatic mesh subdivision technique to extend the scope of use of the proposed method needs to be explored in future work.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare there are no conflicts of interests regarding the publication of this paper.