# An Accelerating Convergence Rate Method for Moving Morphable Components

Authors: Ruichao Lian; Shikai Jing; Zhijun He; Zefang Shi; Guohua Song
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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


---

## Body

## 1. 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–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–12] and level set method [13–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. Then, 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.The 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. Therefore, 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 order to overcome the above difficulties, Guo et al. [17, 18] proposed a so-called moving morphable component method containing precise geometric information, where a series of structural components with explicit geometry information are used as basic building blocks, and the structural topology can be obtained by optimizing the explicit geometry parameters of components (i.e., the lengths, thickness, tilt angle, and center coordinates), such as in Figure 1. Recent years witnessed a growing interest on developing topology optimization methods based on explicit geometry/topology descriptions [19–34].Figure 1
A schematic illustration of the MMC method: (a) the initial layout of the components, (b) the optimization process, and (c) the finally optimized result.
(a)
(b)
(c)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. They 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. Therefore, in order to improve the convergence rate of the MMC method, we propose a hierarchical mesh subdivision method. This 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. This 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. The effectiveness of the proposed method is verified by 2D and 3D numerical examples.


## 2. 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:(1)φwx>0,ifx∈Ωw,φwx=0,ifx∈∂Ωw,φwx<0,ifx∈RΩw,with(2)φwx=maxφ1x,φ2x,…,φix,…,φNx,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(3)φix>0,ifx∈Ωiw,φix=0,ifx∈∂Ωiw,φix<0,ifx∈RΩiw,where Ωiw denotes the design domain comprised by the ith component and Ωw=∪i=1NΩiw, 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].Figure 2
The geometry representation based on the MMC method.According to [17, 18], the geometry of a component is represented by constructing the TDF. The 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. The TDF of the ith component can be denoted as(4)φix,y=a′li6+b′fa′6−1,where(5)fa′=ti1+ti2−2ti32li2a′2+ti2−ti12lia′+ti3,a′b′=cosθisinθi−sinθicosθix−xi0y−yi0,in which xi0, yi0, li, sinθ, ti1, ti2, and ti3 are defined as in Figure 3. For the description of a 3D structural component, see [21].Figure 3
The component with quadratically varying thicknesses.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(6)FindD=D1Τ,D2Τ,…,DiΤ,…,DNΤΤMinimizeJD=UDΤFSubject toKDUD=FVD−Vmax≤0,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, Di,i=1,2,…,N represents the vector of design variables associated with the ith component, where N is the total number of components, and Vmax denotes the upper limit value of VD.


## 3. 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 Figure4, 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. Therefore, the rate of which the topology structure is constructed by component movement is related to the fineness of the mesh.Figure 4
Mesh nodes with distance units: red dots indicate the original mesh nodes and blue open circles indicate mesh nodes with distance units.As shown in Figure5, 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). The 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,ydenotes the center coordinate of one component. Then, the center coordinate position of the adjacent iteration step component can be expressed as(7)xy=xjyjxj+1yj+1,where j represents the iteration step. The moving distances of the components I and J can be calculated by dj=xj+1−xj2+yj+1−yj2. We take j=1,2; the relevant data is shown in Table 1.Figure 5
The moving distance of component I (center is the ring) and J (center is the cross) under different meshes: (a)40×20 FE mesh and (b) 240×120 FE mesh.
(a)
(b)Table 1
Center coordinates and moving distances of components under different meshes.
nel
x
×
nel
y
Component
x
1
,
y
1
x
2
,
y
2
x
3
,
y
3
d
1
d
2
40
×
20
I
(0.75, 0.75)
(0.67, 0.71)
(1.38, 0.29)
0.089
0.825
J
(1.25, 0.25)
(0.60, 0.70)
(1.40, 0.31)
0.791
0.890
240
×
120
I
(0.75, 0.75)
(0.72, 0.74)
(1.29, 1.30)
0.032
0.800
J
(1.25, 0.25)
(0.71, 0.71)
(0.26, 0.29)
0.709
0.616It can be clearly seen from Table1 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. Then, the obtained optimal solution is used to construct a starting point for the next finer mesh, which will speed up the convergence rate.Therefore, a natural idea is to employ the multigrid approaches to the MMC framework. Several multigrid methods have been used in density-based topology optimization. The 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–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)). This 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. Therefore, this paper applies a uniform size hierarchical meshing strategy for the MMC framework.Figure 6
The AMF method: (a) the component that uses the AMF method and (b) movement and display of a component.
(a)
(b)Based on the above analysis, a method of combining a uniform size grid strategy and the MMC framework is proposed. The 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,Pq, that we need to solve can be described as(8)Pq:minfx,s.t.x∈0,1nq,where nq=4nq−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 Pq is solved by the method proposed in this paper, and the obtained optimization solution of Pq is used to construct a starting point for the next subproblem Pq+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=x0,y0,l,t1,t2,t3,sinθ to the variable set of all components: D=D1Τ,D2Τ,…,DiΤ,…,DNΤΤ. 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.Figure 7
The hierarchical mesh refinement strategy: (a) initial layout of components under coarse meshes, (b) optimization results under coarse mesh, (c) topological structure of the coarse mesh results mapped to the finer mesh, and (d) the final topology structure obtained with the finer and finer mesh.
(a)
(b)
(c)
(d)In the hierarchical mesh subdivision method, the number of discrete mesh can be expressed as(9)mesh=nelx×nely,if Iters<Frc,2×nelx×2×nely,if Iters≥Frc,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(10)Frc=Max iternj,where Max iter denotes the maximum number of iterations, nj 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. The number of meshes will be divided into 2j×nelx×2j×nely mesh. As shown in Figure 8, the center coordinates of a component Aik 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. Therefore, after the coarse mesh is subdivided, the equivalent transformation of the topology structure can be realized by mapping the components variables set D′=D1′Τ,D2′Τ,…,Di′Τ,…,DN′ΤΤ with Di′=xi0′,yi0′,li′,ti1′,ti2′,ti3′,sinθi′Τ 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 2j×nelx×2j×nely×2j×nelz, and the same component variables equivalent mapping method can be used to achieve hierarchical mesh subdivision optimization. The biggest difference of the strategy is that the mesh can be subdivided as needed during the optimization process.Figure 8
Equivalent mapping of the 2D component under different meshes.


## 4. Numerical Solution Aspects

### 4.1. Finite Element Analysis

In the paper, the Eulerian mesh is utilized and TDF is used to calculate the value of each node. The 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]:(11)Ee=E0∑j=14Hφjeq4,where q is an integer number, we take 2, E0 denotes Young’s modulus for the material of each element, H=Hx is a regularized Heaviside function, and φje,j=1,…,4, are the values of the TDF function of the whole structure at four nodes of an element:(12)Hx=1,ifx>ε,31−a4×xε−x′3ε′+1+a2,if−ε≤x≤ε,a,ifx<−ε,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(13)ki=Eiks,where ks is the element stiffness matrix when Ei=1.


### 4.2. Sensitivity Analysis

Once the stiffness matrix of each element is been obtained, the global stiffness matrixK can be assembled. The objective function C and the displacement vector U will be constructed as(14)C=UΤKU,F=KU.Assume thatμ is any one of the variables, and its sensitivity of the structural compliance can be written as(15)∂C∂μ=∂UΤKU∂μ=∂UΤ∂μKU+UΤ∂K∂μU+UΤK∂U∂μ,(16)∂F∂μ=∂KU∂μ=∂K∂μU+K∂U∂μ.Since the force is independent of the design variables, we can get (17) by (11) and (13):(17)∂K∂μ=∑i=1N∂∂μki=E4∑i=1N∑j=14qHφjeq−1∂∂μHφjeki.By equations (15), (16), and (17), the derivative of compliance could be written as(18)∂C∂μ=−UΤ∂K∂μU=−UΤE4∑i=1N∑j=14qHφjeq−1∂∂μHφjekiU.For the derivative of volume constraint, we also have(19)∂V∂μ=14∑i=1N∑j=14∂Hφje∂μ.


## 4.1. Finite Element Analysis

In the paper, the Eulerian mesh is utilized and TDF is used to calculate the value of each node. The 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]:(11)Ee=E0∑j=14Hφjeq4,where q is an integer number, we take 2, E0 denotes Young’s modulus for the material of each element, H=Hx is a regularized Heaviside function, and φje,j=1,…,4, are the values of the TDF function of the whole structure at four nodes of an element:(12)Hx=1,ifx>ε,31−a4×xε−x′3ε′+1+a2,if−ε≤x≤ε,a,ifx<−ε,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(13)ki=Eiks,where ks is the element stiffness matrix when Ei=1.


## 4.2. Sensitivity Analysis

Once the stiffness matrix of each element is been obtained, the global stiffness matrixK can be assembled. The objective function C and the displacement vector U will be constructed as(14)C=UΤKU,F=KU.Assume thatμ is any one of the variables, and its sensitivity of the structural compliance can be written as(15)∂C∂μ=∂UΤKU∂μ=∂UΤ∂μKU+UΤ∂K∂μU+UΤK∂U∂μ,(16)∂F∂μ=∂KU∂μ=∂K∂μU+K∂U∂μ.Since the force is independent of the design variables, we can get (17) by (11) and (13):(17)∂K∂μ=∑i=1N∂∂μki=E4∑i=1N∑j=14qHφjeq−1∂∂μHφjeki.By equations (15), (16), and (17), the derivative of compliance could be written as(18)∂C∂μ=−UΤ∂K∂μU=−UΤE4∑i=1N∑j=14qHφjeq−1∂∂μHφjekiU.For the derivative of volume constraint, we also have(19)∂V∂μ=14∑i=1N∑j=14∂Hφje∂μ.


## 5. 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 isE=1 and Poisson’s ratio is υ=0.3. The 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. The 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.


### 5.1. The 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. The design domain, boundary conditions, and external load are shown in Figure9. A unit vertical load is imposed at the middle point of the upper surface. The upper bound of available solid material volume constraint is V¯=0.4. The 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. The initial parameters of each component are all set as 0.380.040.060.040.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. The sizes of the refined meshes that are actually involved in analysis are 120×40 and 240×80 mesh. In (10), j=0,1,2, n0=10, n1=10/3, and n2=5/3. In the SIMP method, the penalty parameter p=3 and the filtering radius r=6.Figure 9
The 2D MBB beam example.Figure 10
The initial component layout of the 2D MBB beam.Some intermediate steps and the history of the convergence process of the three methods are presented in Tables2 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 gray-scale 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. This may not only slow down the convergence of the optimization process, but also affect the final structure. The 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. The topology structure rarely undergoes significant mutations after this solution, which also enhances the stability of the optimization processes at next stages.Table 2
Comparison of topology optimization process using the MMC method and the SIMP approach for the 2D MBB beam.
Approach
MMC
SIMP
Elements
240
×
80
Iterative step 10
Iterative step 50
Iterative step 100
Iterative step 300
—
Convergence historyTable 3
The topology optimization process using the MMC-HMR approach for the 2D MBB beam.
Approach
MMC-HMR
Elements
Level 1:60×20
Level 2:120×40
Level 3:240×80
Iterative step 10
Iterative step 50
Iterative step 100
Iterative step 300
—
Convergence historyThe final optimal results are shown in Table4. 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. The 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. The MMC method does not reach convergence within the set maximum number of iterations, and the total time cost is 2045.26 s. The MMC-HMR method reached convergence at the 639th iterative step. The 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.Table 4
Comparison of the final optimization results of different methods for the 2D MBB beam example.
Method
Elements
Obj
Iters
Time (s)
Optimal topology
MMC
240
×
80
232.38
1000
2045.26
MMC-HMR
228.28
639
709.09
SIMP
258.66
114
17645.88


### 5.2. A 3D Short Cantilever Beam Example

The classical 3D short cantilever beam example is examined. The boundary conditions and external load are shown in Figure11. Totally eighteen components are placed in the design domain. This problem is optimized by the MMC method, the MMC-HMR method, and the SIMP approach, respectively. The 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. The initial components layout is shown in Figure 12. The design domain is discretized by 60×6×30 finite element mesh. In the MMC-HMR method, j=0,1 and n0=5, and the design domain is first discretized into 30×3×15 coarse mesh. When the mesh subdivision condition is met, the 30×3×15 FE mesh will be refined 60×6×30 FE mesh. The upper bound of available solid material volume constraint is V¯=0.3.Figure 11
A 3D short cantilever beam example.Figure 12
The initial design of the 3D short cantilever beam example.The final optimal structural topology obtained by the different methods is shown in Figure13. 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 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. The 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. The intermediate structure topology obtained by the SIMP method also has not changed much, but the disadvantages are similar to the previous 2D optimization result. The structure boundary is nonsmooth and contains more gray elements.Figure 13
The final optimized structure of the 3D short cantilever beam example: (a) the MMC method, (b) the MMC-HMR method, and (c) the SIMP method.
(a)
(b)
(c)Figure 14
Some intermediate designs during the MMC method optimization process.
(a)
(b)Figure 15
Some intermediate designs during the MMC-HMR method optimization process: (a) the optimized structure with a30×3×15 FE mesh and (b) the result of (a) is mapped to the structure with a 60×6×30 FE mesh.
(a)
(b)Figure 16
Some intermediate designs during the SIMP method optimization process.
(a)
(b)Table5 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 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. The 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.Table 5
Comparison of the final optimization results of different methods for the 3D short cantilever example.
Method
Elements
Obj
Iters
Time (s)
Convergence history
MMC
60
×
30
×
6
51.90
4000
42302.46
MMC-HMR
52.22
396
1226.23
SIMP
1292.47
100
2536.48


## 5.1. The 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. The design domain, boundary conditions, and external load are shown in Figure9. A unit vertical load is imposed at the middle point of the upper surface. The upper bound of available solid material volume constraint is V¯=0.4. The 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. The initial parameters of each component are all set as 0.380.040.060.040.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. The sizes of the refined meshes that are actually involved in analysis are 120×40 and 240×80 mesh. In (10), j=0,1,2, n0=10, n1=10/3, and n2=5/3. In the SIMP method, the penalty parameter p=3 and the filtering radius r=6.Figure 9
The 2D MBB beam example.Figure 10
The initial component layout of the 2D MBB beam.Some intermediate steps and the history of the convergence process of the three methods are presented in Tables2 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 gray-scale 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. This may not only slow down the convergence of the optimization process, but also affect the final structure. The 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. The topology structure rarely undergoes significant mutations after this solution, which also enhances the stability of the optimization processes at next stages.Table 2
Comparison of topology optimization process using the MMC method and the SIMP approach for the 2D MBB beam.
Approach
MMC
SIMP
Elements
240
×
80
Iterative step 10
Iterative step 50
Iterative step 100
Iterative step 300
—
Convergence historyTable 3
The topology optimization process using the MMC-HMR approach for the 2D MBB beam.
Approach
MMC-HMR
Elements
Level 1:60×20
Level 2:120×40
Level 3:240×80
Iterative step 10
Iterative step 50
Iterative step 100
Iterative step 300
—
Convergence historyThe final optimal results are shown in Table4. 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. The 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. The MMC method does not reach convergence within the set maximum number of iterations, and the total time cost is 2045.26 s. The MMC-HMR method reached convergence at the 639th iterative step. The 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.Table 4
Comparison of the final optimization results of different methods for the 2D MBB beam example.
Method
Elements
Obj
Iters
Time (s)
Optimal topology
MMC
240
×
80
232.38
1000
2045.26
MMC-HMR
228.28
639
709.09
SIMP
258.66
114
17645.88


## 5.2. A 3D Short Cantilever Beam Example

The classical 3D short cantilever beam example is examined. The boundary conditions and external load are shown in Figure11. Totally eighteen components are placed in the design domain. This problem is optimized by the MMC method, the MMC-HMR method, and the SIMP approach, respectively. The 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. The initial components layout is shown in Figure 12. The design domain is discretized by 60×6×30 finite element mesh. In the MMC-HMR method, j=0,1 and n0=5, and the design domain is first discretized into 30×3×15 coarse mesh. When the mesh subdivision condition is met, the 30×3×15 FE mesh will be refined 60×6×30 FE mesh. The upper bound of available solid material volume constraint is V¯=0.3.Figure 11
A 3D short cantilever beam example.Figure 12
The initial design of the 3D short cantilever beam example.The final optimal structural topology obtained by the different methods is shown in Figure13. 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 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. The 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. The intermediate structure topology obtained by the SIMP method also has not changed much, but the disadvantages are similar to the previous 2D optimization result. The structure boundary is nonsmooth and contains more gray elements.Figure 13
The final optimized structure of the 3D short cantilever beam example: (a) the MMC method, (b) the MMC-HMR method, and (c) the SIMP method.
(a)
(b)
(c)Figure 14
Some intermediate designs during the MMC method optimization process.
(a)
(b)Figure 15
Some intermediate designs during the MMC-HMR method optimization process: (a) the optimized structure with a30×3×15 FE mesh and (b) the result of (a) is mapped to the structure with a 60×6×30 FE mesh.
(a)
(b)Figure 16
Some intermediate designs during the SIMP method optimization process.
(a)
(b)Table5 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 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. The 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.Table 5
Comparison of the final optimization results of different methods for the 3D short cantilever example.
Method
Elements
Obj
Iters
Time (s)
Convergence history
MMC
60
×
30
×
6
51.90
4000
42302.46
MMC-HMR
52.22
396
1226.23
SIMP
1292.47
100
2536.48


## 6. 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. The 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. This 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.



---

*Source: 2478292-2020-01-30.xml*