# Topology Optimization considering Nonsmooth Structural Boundaries in the Intersection Areas of the Components

Authors: Ruichao Lian; Shikai Jing; Zefang Shi; Zhijun He; Guohua Song
Journal: Shock and Vibration
 (2020)

---

## Abstract

In the structural topology optimization approaches, the Moving Morphable Component (MMC) is a new method to obtain the optimized structural topologies by optimizing shapes, sizes, and locations of components. However, the optimized structure boundary usually generates local nonsmooth areas due to incomplete connection between components. In the present paper, a topology optimization approach considering nonsmooth structural boundaries in the intersection areas of the components based on the MMC is proposed. The variability of components’ shape can be obtained by constructing the topology description function (TDF) with multiple thickness and length variables. The shape of components can be modified according to the structural responses during the optimization process, and the relatively smooth structural boundaries are generated in the intersection areas of the components. To reduce the impact of the initial layout on the rate of convergence, this method is implemented in a hierarchical variable calling strategy. Compared with the original MMC method, the advantage of the proposed approach is that the smoothness of the structural boundaries can be effectively improved and the geometric modeling ability can be enhanced in a concise way. The effectiveness of the proposed method is demonstrated for topology optimization of the minimum compliance problem and compliant mechanisms.


---

## Body

## 1. Introduction

Structural topology optimization aims to find the best distribution of materials within a prescribed design domain using an optimization algorithm in order to achieve some exceptional structural performance [1]. It breaks the limitations of the designer’s thinking in traditional structure design and can get a more liberal, updated, and lighter conceptual design. At the same time, structural topology optimization is also considered a more challenging design approach compared with the size optimization and shape optimization. Since Bendsøe and Kikuchi [2] developed the topology optimization method based on the homogenization theory, after decades of development and improvement, topology optimization has gradually become an important method in the structural concept design stage and has been applied successfully to a wide range of physical disciplines [3–6]. So far, the solid isotropic material with penalization (SIMP) method [7–9] and level set method [10–14] are two popular strategies. For the other methods, the reader is referred to [15–17]. In the SIMP approach, the design domain is discretized by finite element mesh. Then the pseudodensity of each element is set as a design variable to characterize the form of material distribution. It transforms the structural topology optimization problems into material distribution problems. The main advantages include simpler concepts, high computational efficiency, and easy to implement [18–20]. However, the densities of a certain number of elements are left to represent the final optimal results; the structure boundaries are easy to generate with many grey-scale elements and jagged shapes [21–23]. In the level set method, the values of the level set function at the node points in the finite element mesh are often taken as the design variables. The structural boundaries can be specified by drawing the contour of the level set function in one-higher dimensional space. The optimization process is represented by the evolution of specific faces. Therefore, it often obtains the optimization result with smooth boundaries and easily handles topology changes compared to the traditional topology optimization methods [24–26]. But whether the SIMP approach or level set method, the design variables increase sharply with the mesh refinement, especially in the three-dimension problem, which will far exceed the solution capacity of the existing optimization algorithm [27]. In SIMP and level set methods, an implicit way is used to solve topology optimization problems; that is, without explicitly geometric information, it is implanted. It is difficult to accurately control the structural feature sizes, which is an important factor in considering manufacturability [28]. At the same time, the geometric/topology representation is not quite consistent with that in modern computer-aided design/engineering (CAD/CAE) modeling systems [29, 30]. This treatment sets up a barrier in optimizing structural dimensions and direct link between the final topology results and CAD/CAE systems [31].In order to overcome the above problems, recently, Guo et al. [27] proposed a more accurate and geometric topology optimization method based on the topology optimization framework, a so-called Moving Morphable Component (MMC). In this framework, a series of components with movable and deformable capabilities are used as building blocks of the topology optimization, and the optimal structural topology that meets specific performance can be obtained by varying the lengths, thickness, tilt angle, and center coordinates of these components. The same idea has also been applied in continuum-based topology optimization with discrete elements [32–36]. Although any curved structural parts can be approximated by precisely controlling a certain number of components with uniform thicknesses [27], so as to enhance the smoothness of the structural boundaries and the geometry modeling capability of the method, it is always highly desirable to include components with variable shapes in the MMC approach. Recently, Zhang et al. [37, 38] improved the geometry modeling capability of the MMC approach by constructing the TDF with variable thickness components. Guo et al. [39] proposed the TDF of the component with a curved skeleton to constructed curved structural parts with a smooth boundary. In Meisam Takalloozadeh et al.’s method [40], the start and end of the component with two straight lines were modified as a concave/convex pointed shape by adding a new variable. Van-Nam Hoang et al. [33] and Deng et al. [32] constructed the rectangle moving component with semicircular ends and enhanced connectivity between components by adding constraints.Although the geometric modeling capabilities of the above methods and the connection methods of components have been studied, in [37, 40], there are still some obvious nonsmooth boundaries in the optimized structure using the straight lines or concave/convex pointed shapes at the start and end of a component. The method of [39] finds the optimal structural by constructing the components with curved skeleton and most of the boundaries of the optimized structure are smooth, but just like as described at the end of their paper, there still exist some local nonsmooth boundaries especially in the intersection areas of the components. In the optimization process of [32, 33], the thickness of the component changed uniformly, which limits the geometry modeling capability of their methods.Therefore, based on the MMC framework, this paper mainly studies how to use the variability of components’ shape to solve the local nonsmooth problem in the intersection areas of the components. For this purpose, we control the length of components by adding new design variables in the TDF of the quadratically varying thicknesses. The variability of components’ shape and the geometry modeling capability of the original MMC method are enhanced through the coordination of length and thickness variables. In the topology optimization process, the shape of components can be changed according to the structural responses, and relatively smooth structural boundaries are generated in the intersection areas between the components. In addition, considering the influence of the initial layout with more components on the rate of the convergence, the hierarchical variables calling strategy is proposed to optimize the layout of components in the initial stage. The effectiveness of the proposed method is verified by several numerical examples.


## 2. Problem Formulation

In the MMC approach, the topology description function (TDF,φx) of components can be expressed in the following form:(1)φx>0,ifx∈Ω,φx=0,ifx∈∂Ω,φx<0,ifx∈R\Ω,where x represents a set of points within the design domain; R represents the entire design area; Ω⊂R denotes a subset of R occupied by N components made of the physical material. If the number of components included in the design area is N, the TDF of the whole solid material in the structure area can be built as(2)φx=maxφ1x,…,φix,…,φNx,where i=1,…,N is the number of components and the TDF of the i−th component can be written as(3)φix>0,ifx∈Ωi,φix=0,ifx∈∂Ωi,φix<0,ifx∈R\Ωi,where Ωi is the area occupied by the i−th and Ω=∪i=1NΩi. Then the geometric representation of the corresponding structure is shown in Figure 1.Figure 1
The geometry representation based on the MMC framework.Based on the MMC topology optimization framework, a typical topology optimization problem can be described as(4)Find,R=R1,…,Ri,…,RN,to minimize,J=JR,subject to,ξkR≤0,k=1,…,m,R⊂UR,where JR and ξk, k=1,…,m, denote the objective function and constraint function, respectively, and UR is the admissible set that R belongs to.Topology optimization for a mean compliance minimization under available volume constraint problem is considered, and the corresponding mathematical model can be written as(5)find,R=R1,…,Ri,…,RN,ux∈GΩS,to minimize,J=∫RHφSx;Rf⋅udV+∫Γtt⋅udS,subject to,∫RHqφSx;RE:εu:ευdV=∫RHφSx;Rf⋅υdV+∫Γtt⋅υdS,∀υ∈UI,∫RHφSx;RdV≤V¯,R⊂UR,u=u¯,onΓu,where Ri, i=1,…,N, is the vector of the design variables of the i−th component. J is the objective function. f, u, t, and ε are the body force density, the displacement field, the prescribed surface traction on Neumann boundary Γt, and the linear strain tensor, respectively. u¯ is the displacement on Dirichlet boundary Γu. υ is the corresponding test function defined on UI=υ|υ∈HΩS,υ=0onΓu, and the symbol H=Hx means the Heaviside function. φSx;R denotes the topology description function set for the overall component. q>1 is an integer (in this work, q=2 is used), and E=ES/1+vSΠ+vS/1−2vSδ⊗δ is the isotropic elastic tensor. E and v are Young’s modulus and Poisson’s ratio of the solid material, respectively, and Π and δ denote the symmetric part of the fourth-order identity tensor and the second-order identity tensor, respectively. V¯ is the upper limit of the volume of the solid material.


## 3. TDF for the Variability of Components’ Shape and Hierarchical Variable Calling Strategy

In this part, we will discuss how to construct the TDF for the variability of components’ shape and implement the hierarchical variables calling strategy based on the MMC framework.


### 3.1. TDF for the Variability of Components’ Shape

The difference consists that, in the MMC method, it is possible to give an explicit description of the boundary and geometry features of a component. This cannot be achieved in the conventional level set method. In [27], the uniform thickness component can be described using the hyperelliptic TDF. In order to improve the geometry modeling capabilities of the MMC approach, the TDF of quadraticall varying thicknesses (QVT) can be expressed as(6)φix,y=x′lm+y′fx′m−1,where m is a relatively large positive integer, which takes (6) in this paper. Here, fx′ controls the shape of components(7)fx′=ti1+ti2−2ti32li2x′2+ti2−ti12lix′+ti3,(8)x′y′=cosθisinθi−sinθicosθix−x0iy−y0i,where x0i,y0i, li, and θi denote the coordinate of the center, the half-length, and the inclined angle of i−th component, respectively. ti1, ti2, and ti3 are the three thicknesses for different positions of the i−th component, respectively, as shown in Figure 2.Figure 2
Geometry description of the quadratically variable component.It is worth noting that, in this method, the start and end of the component were two straight lines perpendicular to the direction of the component. When two or more components are connected at an angle, as shown in Figure3 which plots the components with dashed boundary, the local nonsmooth region of the structural boundary is often caused. If we assume that the straight line of the start and end of the component can be bent and deformed, the variability of the shape of the component can be achieved by combining the thickness variations. Then, the local smooth structural boundaries can be formed in the intersection parts of the components according to the structural response in the optimization process, as shown in Figure 3 which plots the components with solid boundary. The geometry description of the components’ shape with variability is shown in Figure 4. The deformation of components’ shape can be achieved by adjusting the length and thickness variables. The corresponding TDF can be constructed as(9)ψix,y=x′gy′m+y′fx′m−1,where fx′ is the same as equation (7) and gy′ can be constructed as(10)gy′=r2−y′−y0is2+x0is,ifli3<li1+li22,li2+li12+li2−li12ti1y′,ifli3=li1+li22,−r2−y′−y0is2+x0is,ifli3>li1+li22,where li1, li2, and li3 are the three lengths for different positions of a component as shown in Figure 4, respectively. x0is,y0is and r are the center coordinates and radius of the arc of the start and end of a component, respectively, and we have(11)x0is=li12+ti12−li22+ti12−2li322li3−li1+li3−li2,y0is=li3−li2li12+ti122li3−li1+li3−li2+li3−li1li22+ti12+li3li1−li22li3−li1+li3−li2,r2=x0is+li12+y0is+ti12.Figure 3
The smooth boundaries construction.Figure 4
Geometry description of the variability of components’ shape.By inputting different parameters in equation (9), the variability of components’ shape can be easily achieved. Some components are shown in Figure 5, and the values for the lengths and thicknesses variables l1l2l3t1t2t3 of components are taken as (a) [0.6 0.6 0.9 0.3 0.3 0.3], (b) [1.0 1.0 0.7 0.3 0.3 0.3], (c) [1.0 0.4 0.7 0.3 0.3 0.3], (d) [1.0 0.1 0.8 0.3 0.3 0.3], (e) [1.0 0.05 0.7 0.08 0.5 0.15], and (f) [1.0 0.1 0.7 0.7 0.7 0.03], respectively. It is worth noting that although the topology optimization process using the variability of components’ shape based on the MMC framework is somewhat similar to the level set method (TDF has been used to represent the geometry of components), the proposed method can also possibly give an explicit description of the boundary and geometric features (such as length and thickness) of components, which cannot be achieved in the level set method that uses the free transformation of the level set functions to represent the boundary of a structure.Figure 5
Geometry description of components with different parameters.
(a)(b)(c)(d)(e)(f)


### 3.2. Hierarchical Variables Calling Strategy

In this section, we will discuss how to optimize the initial layout of components to improve the rate of convergence. In MMC framework, although structural topology can be constructed by moving, deforming, disappearing, and overlapping of a set of components, the different initial layouts of components have different effects on the rate of convergence. In particular, when there are more components in the design domain, the effective force transfer path may be affected by improper positioning or deformation of some components, which will delay the rate of convergence.In order to improve the rate of convergence of the optimization process, a hierarchical variables calling strategy is proposed. In this strategy, the optimization process is divided into two levels: the level 1 is the layout variables calling and the level 2 is all variables calling. The average variation of components’ positioning (the central coordinates and the rotation angles) within the adjacent iterations is set as a constraint. In the initial stage of optimization, only the layout variables (the central coordinates and the rotation angles) are called, and the initial layout of the components can be modified according to the structural responses in this process. When the constraint is satisfied, the full variables are called to construct the topology optimization structure. It is worth noting that the output of the level 1 is used as input for the level 2.The constraint of the hierarchical variables calling strategy can be expressed as(12)Sd_acj+1−Sd_acj/Sd_acj−Sd_acj−1≤τ,Sζ_acj+1−Sζ_acj/Sζ_acj−Sζ_acj−1≤σ,where j is the j−th iteration, τ and σ are the set maximum upper limits of change values, respectively, and generally, the smaller positive numbers are selected. Sdacj and Sζacj are the average change values of the center coordinates and rotation angle of components:(13)Sdacj=∑i=1NdtopiN,Sζacj=∑i=1NζtcpiN,with(14)dtopi=x0ij+1−x0ij2+y0ij+1−y0ij2,ζtcpi=ζij+1−ζij,where dtopi and ζtcpi are the center coordinates x0i,y0i and the rotation angles ζi change value of the i−th component in twice adjacent iterations, respectively, as shown in Figure 6.Figure 6
The center coordinates and angle changes.


## 3.1. TDF for the Variability of Components’ Shape

The difference consists that, in the MMC method, it is possible to give an explicit description of the boundary and geometry features of a component. This cannot be achieved in the conventional level set method. In [27], the uniform thickness component can be described using the hyperelliptic TDF. In order to improve the geometry modeling capabilities of the MMC approach, the TDF of quadraticall varying thicknesses (QVT) can be expressed as(6)φix,y=x′lm+y′fx′m−1,where m is a relatively large positive integer, which takes (6) in this paper. Here, fx′ controls the shape of components(7)fx′=ti1+ti2−2ti32li2x′2+ti2−ti12lix′+ti3,(8)x′y′=cosθisinθi−sinθicosθix−x0iy−y0i,where x0i,y0i, li, and θi denote the coordinate of the center, the half-length, and the inclined angle of i−th component, respectively. ti1, ti2, and ti3 are the three thicknesses for different positions of the i−th component, respectively, as shown in Figure 2.Figure 2
Geometry description of the quadratically variable component.It is worth noting that, in this method, the start and end of the component were two straight lines perpendicular to the direction of the component. When two or more components are connected at an angle, as shown in Figure3 which plots the components with dashed boundary, the local nonsmooth region of the structural boundary is often caused. If we assume that the straight line of the start and end of the component can be bent and deformed, the variability of the shape of the component can be achieved by combining the thickness variations. Then, the local smooth structural boundaries can be formed in the intersection parts of the components according to the structural response in the optimization process, as shown in Figure 3 which plots the components with solid boundary. The geometry description of the components’ shape with variability is shown in Figure 4. The deformation of components’ shape can be achieved by adjusting the length and thickness variables. The corresponding TDF can be constructed as(9)ψix,y=x′gy′m+y′fx′m−1,where fx′ is the same as equation (7) and gy′ can be constructed as(10)gy′=r2−y′−y0is2+x0is,ifli3<li1+li22,li2+li12+li2−li12ti1y′,ifli3=li1+li22,−r2−y′−y0is2+x0is,ifli3>li1+li22,where li1, li2, and li3 are the three lengths for different positions of a component as shown in Figure 4, respectively. x0is,y0is and r are the center coordinates and radius of the arc of the start and end of a component, respectively, and we have(11)x0is=li12+ti12−li22+ti12−2li322li3−li1+li3−li2,y0is=li3−li2li12+ti122li3−li1+li3−li2+li3−li1li22+ti12+li3li1−li22li3−li1+li3−li2,r2=x0is+li12+y0is+ti12.Figure 3
The smooth boundaries construction.Figure 4
Geometry description of the variability of components’ shape.By inputting different parameters in equation (9), the variability of components’ shape can be easily achieved. Some components are shown in Figure 5, and the values for the lengths and thicknesses variables l1l2l3t1t2t3 of components are taken as (a) [0.6 0.6 0.9 0.3 0.3 0.3], (b) [1.0 1.0 0.7 0.3 0.3 0.3], (c) [1.0 0.4 0.7 0.3 0.3 0.3], (d) [1.0 0.1 0.8 0.3 0.3 0.3], (e) [1.0 0.05 0.7 0.08 0.5 0.15], and (f) [1.0 0.1 0.7 0.7 0.7 0.03], respectively. It is worth noting that although the topology optimization process using the variability of components’ shape based on the MMC framework is somewhat similar to the level set method (TDF has been used to represent the geometry of components), the proposed method can also possibly give an explicit description of the boundary and geometric features (such as length and thickness) of components, which cannot be achieved in the level set method that uses the free transformation of the level set functions to represent the boundary of a structure.Figure 5
Geometry description of components with different parameters.
(a)(b)(c)(d)(e)(f)


## 3.2. Hierarchical Variables Calling Strategy

In this section, we will discuss how to optimize the initial layout of components to improve the rate of convergence. In MMC framework, although structural topology can be constructed by moving, deforming, disappearing, and overlapping of a set of components, the different initial layouts of components have different effects on the rate of convergence. In particular, when there are more components in the design domain, the effective force transfer path may be affected by improper positioning or deformation of some components, which will delay the rate of convergence.In order to improve the rate of convergence of the optimization process, a hierarchical variables calling strategy is proposed. In this strategy, the optimization process is divided into two levels: the level 1 is the layout variables calling and the level 2 is all variables calling. The average variation of components’ positioning (the central coordinates and the rotation angles) within the adjacent iterations is set as a constraint. In the initial stage of optimization, only the layout variables (the central coordinates and the rotation angles) are called, and the initial layout of the components can be modified according to the structural responses in this process. When the constraint is satisfied, the full variables are called to construct the topology optimization structure. It is worth noting that the output of the level 1 is used as input for the level 2.The constraint of the hierarchical variables calling strategy can be expressed as(12)Sd_acj+1−Sd_acj/Sd_acj−Sd_acj−1≤τ,Sζ_acj+1−Sζ_acj/Sζ_acj−Sζ_acj−1≤σ,where j is the j−th iteration, τ and σ are the set maximum upper limits of change values, respectively, and generally, the smaller positive numbers are selected. Sdacj and Sζacj are the average change values of the center coordinates and rotation angle of components:(13)Sdacj=∑i=1NdtopiN,Sζacj=∑i=1NζtcpiN,with(14)dtopi=x0ij+1−x0ij2+y0ij+1−y0ij2,ζtcpi=ζij+1−ζij,where dtopi and ζtcpi are the center coordinates x0i,y0i and the rotation angles ζi change value of the i−th component in twice adjacent iterations, respectively, as shown in Figure 6.Figure 6
The center coordinates and angle changes.


## 4. Numerical Solution Aspects

In this section, numerical solution aspects in the proposed approach will be discussed in detail.


### 4.1. Finite Element Analysis

In this study, in order to analyze the structure constructed using the MMC method, the Eulerian mesh is utilized and TFD is used to calculate each node. In order to enhance the computational efficiency, the ersatz material model is used for FEM analysis. Young’s modulus of every element can be calculated by equation (15), which was proposed in [37]:(15)Ee=E0∑j=14Hφjeq4,where q is an integer number, we take q=2, E0 is Young’s modulus for the material of every element, and H is a regularized Heaviside function:(16)Hx=1,ifε<x,31−a4×xε−x′3ε′+1+a2,otherwise,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 i−th element can be obtained by equation (17):(17)ki=EikS.In the above formula,kS is the stiffness matrix of the element, regardless of the number of materials occupied by components.


### 4.2. Sensitivity Analysis

Only design sensitivity analysis for the compliance minimization problem is discussed here. Other objectives/constraints can be achieved by resorting to adjoint sensitivity analysis. When we know the stiffness matrix of each element, the global stiffness matrix can be assembled. The structural compliance and the displacement vector will be obtained by using equation (18):(18)C=UTKU,F=KU.Whenμ is assumed to be any one of the variables, then the sensitivity of the structural compliance of μ can be expressed as(19)∂C∂μ=∂UTKU∂μ=∂UT∂μKU+UT∂K∂μU+UTK∂U∂μ,(20)∂F∂μ=∂KU∂μ=∂K∂μU+K∂U∂μ.Since the force is independent of the design variables, equation (20) can be transformed into(21)∂K∂μU=−K∂U∂μ.By equations (15) and (17), we can get equation (22):(22)∂K∂μ=∑i=1N∂∂μki=E4∑i=1N∑j=14qHφjeq−1∂∂μHφjeki.By substituting equations (21) and (22) into equation (19), which can be written as follows:(23)∂C∂μ=−UT∂K∂μU=−UTE4∑i=1N∑j=14qHφjeq−1∂∂μHφjekiU,where ki is the element stiffness matrix corresponding to φje=1,j=1,…,4. Because of φix is a function of μ, the derivative of the structural compliance can be calculated using equation (23).For the derivative of volume constraint, we also have(24)∂V∂μ=14∑i=1N∑j=14∂Hφje∂μ,where ∂Hφje/∂μ can be calculated easily since φx is an explicit function of μ.It is worth noting that, in the above derivations, as mentioned in [37], the nondifferentiable issue arising from the max operation when φ is constructed by multiple components overlapping is neglected. The numerical case also proves that this has no effect on the optimization process.


## 4.1. Finite Element Analysis

In this study, in order to analyze the structure constructed using the MMC method, the Eulerian mesh is utilized and TFD is used to calculate each node. In order to enhance the computational efficiency, the ersatz material model is used for FEM analysis. Young’s modulus of every element can be calculated by equation (15), which was proposed in [37]:(15)Ee=E0∑j=14Hφjeq4,where q is an integer number, we take q=2, E0 is Young’s modulus for the material of every element, and H is a regularized Heaviside function:(16)Hx=1,ifε<x,31−a4×xε−x′3ε′+1+a2,otherwise,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 i−th element can be obtained by equation (17):(17)ki=EikS.In the above formula,kS is the stiffness matrix of the element, regardless of the number of materials occupied by components.


## 4.2. Sensitivity Analysis

Only design sensitivity analysis for the compliance minimization problem is discussed here. Other objectives/constraints can be achieved by resorting to adjoint sensitivity analysis. When we know the stiffness matrix of each element, the global stiffness matrix can be assembled. The structural compliance and the displacement vector will be obtained by using equation (18):(18)C=UTKU,F=KU.Whenμ is assumed to be any one of the variables, then the sensitivity of the structural compliance of μ can be expressed as(19)∂C∂μ=∂UTKU∂μ=∂UT∂μKU+UT∂K∂μU+UTK∂U∂μ,(20)∂F∂μ=∂KU∂μ=∂K∂μU+K∂U∂μ.Since the force is independent of the design variables, equation (20) can be transformed into(21)∂K∂μU=−K∂U∂μ.By equations (15) and (17), we can get equation (22):(22)∂K∂μ=∑i=1N∂∂μki=E4∑i=1N∑j=14qHφjeq−1∂∂μHφjeki.By substituting equations (21) and (22) into equation (19), which can be written as follows:(23)∂C∂μ=−UT∂K∂μU=−UTE4∑i=1N∑j=14qHφjeq−1∂∂μHφjekiU,where ki is the element stiffness matrix corresponding to φje=1,j=1,…,4. Because of φix is a function of μ, the derivative of the structural compliance can be calculated using equation (23).For the derivative of volume constraint, we also have(24)∂V∂μ=14∑i=1N∑j=14∂Hφje∂μ,where ∂Hφje/∂μ can be calculated easily since φx is an explicit function of μ.It is worth noting that, in the above derivations, as mentioned in [37], the nondifferentiable issue arising from the max operation when φ is constructed by multiple components overlapping is neglected. The numerical case also proves that this has no effect on the optimization process.


## 5. Numerical Examples

In this section, the effectiveness of the proposed method is verified by several representative numerical examples. Unless otherwise stated, all examples use a four-node bilinear square unit discrete design area. Young’s modulus and Poisson’s ratio of the solid material are taken asE=1 and υ=0.3, respectively. The two-dimensional topology optimization problem is considered. The Method of Moving Asymptotes (MMA) is used in this method, and the hyperparameters of MMA are epsimin=10−10, raa0=0.01, albefa=0.4, asyinit=0.1, asyincr=0.8, and asydecr=0.6. In addition, all the example codes analysis is performed on the LENOVO-TY5060 workstation with an Intel® CoreTm i5-6400 CPU 2.70 GHz, 4G RAM, and Windows10 and ran with MATLAB R2017a.


### 5.1. The Compliant Mechanism Problem

To demonstrate the effectiveness of the proposed method in dealing with local nonsmoothness problems in the intersection areas of the components, the compliant mechanism design problem is considered. The design objective is to minimize the geometric advantage of the compliant mechanism. The design domain, geometry data, and boundary conditions are shown in Figure7. In this problem, the objective is to minimize the geometric advantage (GA) of the compliant mechanism GA=−Δout/Δin, where Δin=Finu1in+Foutu2in and Δout=Finu1out+Foutu2out, ujin,ujout,j=1,2 meaning the displacement in the horizontal direction under a unit force applied at point j. The elastic constant of the springs is kout=0.1. The design domain with length and width of 2×2 is discretized by a 80×80 FE mesh. The upper bound of volume constraint is set to 0.3. Since the problem under consideration is symmetric in nature, only half of the structure is used to optimize in this example. Sixteen components are placed in the initial design domain, and the overall number of design variables is 16×9=144. The initial values of the design variables for each component of the QVT method and this paper method are set as 0.380.040.060.040.7 and 0.390.390.390.040.060.040.7, respectively. The initial layout of the components is shown in Figure 8.Figure 7
The compliant mechanism problem.Figure 8
The initial designs for the compliant mechanism problem: (a) QVT and (b) the paper method.
(a)(b)The intermediate iteration steps and the final topology optimization results for the two methods are shown in Figures9 and 10. It can be seen that the final structure obtained by the QVT method forms an obvious incomplete connection phenomenon in the intersection area between the components, which easily cause the local structural boundaries nonsmoothness (see Figure 9(d)). From the contour plot of the final structures shown in Figure 10, the components exhibit some flexibility in the form of the connection, and a relatively smooth structural boundary transition is obtained. In the case of the proposed method, the hinge-like connection can no longer be produced in the final optimized structure as illustrated in Figure 10(d).Figure 9
Some intermediate optimization results and the final results of the compliant mechanism example with the QVT method.
(a)(b)(c)(d)Figure 10
Some intermediate optimization results and the final results of the compliant mechanism example with the paper method.
(a)(b)(c)(d)Figure11 provides the convergence history of the objective function and the constraint function. It can be observed from Figure 11(b) that, for this example, the value of the objective function has dropped sharply in the earlier iteration step, which indicates that an effective loading transmission path can be faster constructed using the method. The comparison results are shown in Table 1. The values of the objective functions obtained by the two methods are Obj=−0.811 and Obj=−0.851, respectively. Compared with the QVT method, the proposed method has a faster convergence rate and less time consumption.Figure 11
Convergence history of the compliant mechanism example: (a) the QVT method and (b) the paper method.
(a)(b)Table 1
The data obtained from the compliant mechanism problem.
MethodIterationsObjections valueTime (s)The QVT method1000–0.811432.65The paper method716–0.851386.38


### 5.2. The MMB Problem

In the literature on structural optimization studies, this example is often used to verify the validity of a topology optimization method. The design domain and boundary conditions are shown in Figure12. The unit vertical load is imposed on the middle point of the upper surface of the beam.Figure 12
The MBB problem.The design domain with length and width of6×1 is discretized by a 240×40 FE mesh. The upper bound of volume constraint is set to 0.4. Because of the symmetric property of the problem, only half of the structure is used to be optimized in this example.Twenty-four components are placed in the initial design domain, and the overall number of design variables is only24×9=216. The initial values of the design variables for each component of the QVT method and this paper method are set as 0.380.040.060.040.7 and 0.390.410.390.040.060.040.7, respectively. In equation (12), the maximum upper limits of change values are τ=0.38 and σ=0.57. The initial layout of the components is shown in Figure 13.Figure 13
The initial designs for the MBB problem: (a) the QVT method and (b) the paper method.
(a)(b)Some steps of the optimization and the convergence histories of the MBB example are shown in Figures14–16, respectively. From Figure 14(f) and the 810th step in Figure 16, we can observe that the final topology optimization structures obtained by the two methods are very similar.Figure 14
Some intermediate optimization results and the final results of the compliant mechanism example with the QVT method.
(a)(b)(c)(d)(e)(f)Figure 15
Convergence history of the MBB example with the QVT method.Figure 16
Convergence history of the MBB example using the paper method.However, it is worth noting that the QVT method shows a strong instability in the optimization process. In contrast, the proposed method has a more stable convergence process. It can be clearly seen that when equation (12) is not satisfied, the shape and size of the components do not change, as shown the 20th step contour plot in Figure 16. The entire process only optimizes the center coordinates and the rotation angles of components before the 40th iteration. Once equation (12) is satisfied, the components can be moved and deformed according to the structural response, as shown in the contour plot of the 41th and 42th step in Figure 16.The comparison results by using the two methods are presented in Table2. The values of the objective function are Obj=232.84 and Obj=228.59, respectively. Compared with the QVT method, the value of the objective function obtained is slightly lower and the time consumption has been reduced with this paper method. The optimization process has been converged at the 810th iteration. It is worth mentioning that the proposed method enhances the stability of the solution process and accelerates the convergence rate.Table 2
The data obtained from the MBB problem.
MethodIterationsObjections valueTime (s)The QVT method1000232.84584.13The paper method810228.49566.68


### 5.3. The Bridge Problem

In order to demonstrate that this paper method cannot only solve the structural boundaries nonsmoothness in the intersection areas of components but also realize the variability of components’ shape to some extent, the bridge problem is considered in this section. The loading and boundary conditions are shown in Figure17. The unit vertical load is evenly imposed on the upper surface of the design field. The design domain with length and width of 3×1 is discretized by a 150×50 FE mesh. The upper bound of volume constraint is set to 0.4. A rectangle domain with length and width of 3×0.1 on the top of the design domain is a nondesign area.Figure 17
The bridge problem.The initial layout is shown in Figure18. The initial values of the design variables for each component of the QVT method and this paper method are set as 0.380.040.060.040.7 and 0.390.410.390.040.060.040.7, respectively.Figure 18
The initial designs for the bridge problem: (a) the QVT method and (b) the paper method.
(a)(b)From the intermediate steps of the optimization process shown in Figures19 and 20, the components of this paper method have more flexible deformation ability in the optimization process, and the smooth optimization structural boundaries with circular arc degree are obtained. Although the number of components in the design domain is small, the effective load transmission paths can still be obtained by proper bending deformation of the components. The objective function values are Obj=155249.20 and Obj=146920.09, respectively. It can be observed that the curved bridge structures have a smaller objective function value.Figure 19
Some intermediate optimization results and the final results of the bridge example with the QVT method.
(a)(b)(c)(d)Figure 20
Some intermediate optimization results and the final results of the bridge example with the paper method.
(a)(b)(c)(d)


## 5.1. The Compliant Mechanism Problem

To demonstrate the effectiveness of the proposed method in dealing with local nonsmoothness problems in the intersection areas of the components, the compliant mechanism design problem is considered. The design objective is to minimize the geometric advantage of the compliant mechanism. The design domain, geometry data, and boundary conditions are shown in Figure7. In this problem, the objective is to minimize the geometric advantage (GA) of the compliant mechanism GA=−Δout/Δin, where Δin=Finu1in+Foutu2in and Δout=Finu1out+Foutu2out, ujin,ujout,j=1,2 meaning the displacement in the horizontal direction under a unit force applied at point j. The elastic constant of the springs is kout=0.1. The design domain with length and width of 2×2 is discretized by a 80×80 FE mesh. The upper bound of volume constraint is set to 0.3. Since the problem under consideration is symmetric in nature, only half of the structure is used to optimize in this example. Sixteen components are placed in the initial design domain, and the overall number of design variables is 16×9=144. The initial values of the design variables for each component of the QVT method and this paper method are set as 0.380.040.060.040.7 and 0.390.390.390.040.060.040.7, respectively. The initial layout of the components is shown in Figure 8.Figure 7
The compliant mechanism problem.Figure 8
The initial designs for the compliant mechanism problem: (a) QVT and (b) the paper method.
(a)(b)The intermediate iteration steps and the final topology optimization results for the two methods are shown in Figures9 and 10. It can be seen that the final structure obtained by the QVT method forms an obvious incomplete connection phenomenon in the intersection area between the components, which easily cause the local structural boundaries nonsmoothness (see Figure 9(d)). From the contour plot of the final structures shown in Figure 10, the components exhibit some flexibility in the form of the connection, and a relatively smooth structural boundary transition is obtained. In the case of the proposed method, the hinge-like connection can no longer be produced in the final optimized structure as illustrated in Figure 10(d).Figure 9
Some intermediate optimization results and the final results of the compliant mechanism example with the QVT method.
(a)(b)(c)(d)Figure 10
Some intermediate optimization results and the final results of the compliant mechanism example with the paper method.
(a)(b)(c)(d)Figure11 provides the convergence history of the objective function and the constraint function. It can be observed from Figure 11(b) that, for this example, the value of the objective function has dropped sharply in the earlier iteration step, which indicates that an effective loading transmission path can be faster constructed using the method. The comparison results are shown in Table 1. The values of the objective functions obtained by the two methods are Obj=−0.811 and Obj=−0.851, respectively. Compared with the QVT method, the proposed method has a faster convergence rate and less time consumption.Figure 11
Convergence history of the compliant mechanism example: (a) the QVT method and (b) the paper method.
(a)(b)Table 1
The data obtained from the compliant mechanism problem.
MethodIterationsObjections valueTime (s)The QVT method1000–0.811432.65The paper method716–0.851386.38


## 5.2. The MMB Problem

In the literature on structural optimization studies, this example is often used to verify the validity of a topology optimization method. The design domain and boundary conditions are shown in Figure12. The unit vertical load is imposed on the middle point of the upper surface of the beam.Figure 12
The MBB problem.The design domain with length and width of6×1 is discretized by a 240×40 FE mesh. The upper bound of volume constraint is set to 0.4. Because of the symmetric property of the problem, only half of the structure is used to be optimized in this example.Twenty-four components are placed in the initial design domain, and the overall number of design variables is only24×9=216. The initial values of the design variables for each component of the QVT method and this paper method are set as 0.380.040.060.040.7 and 0.390.410.390.040.060.040.7, respectively. In equation (12), the maximum upper limits of change values are τ=0.38 and σ=0.57. The initial layout of the components is shown in Figure 13.Figure 13
The initial designs for the MBB problem: (a) the QVT method and (b) the paper method.
(a)(b)Some steps of the optimization and the convergence histories of the MBB example are shown in Figures14–16, respectively. From Figure 14(f) and the 810th step in Figure 16, we can observe that the final topology optimization structures obtained by the two methods are very similar.Figure 14
Some intermediate optimization results and the final results of the compliant mechanism example with the QVT method.
(a)(b)(c)(d)(e)(f)Figure 15
Convergence history of the MBB example with the QVT method.Figure 16
Convergence history of the MBB example using the paper method.However, it is worth noting that the QVT method shows a strong instability in the optimization process. In contrast, the proposed method has a more stable convergence process. It can be clearly seen that when equation (12) is not satisfied, the shape and size of the components do not change, as shown the 20th step contour plot in Figure 16. The entire process only optimizes the center coordinates and the rotation angles of components before the 40th iteration. Once equation (12) is satisfied, the components can be moved and deformed according to the structural response, as shown in the contour plot of the 41th and 42th step in Figure 16.The comparison results by using the two methods are presented in Table2. The values of the objective function are Obj=232.84 and Obj=228.59, respectively. Compared with the QVT method, the value of the objective function obtained is slightly lower and the time consumption has been reduced with this paper method. The optimization process has been converged at the 810th iteration. It is worth mentioning that the proposed method enhances the stability of the solution process and accelerates the convergence rate.Table 2
The data obtained from the MBB problem.
MethodIterationsObjections valueTime (s)The QVT method1000232.84584.13The paper method810228.49566.68


## 5.3. The Bridge Problem

In order to demonstrate that this paper method cannot only solve the structural boundaries nonsmoothness in the intersection areas of components but also realize the variability of components’ shape to some extent, the bridge problem is considered in this section. The loading and boundary conditions are shown in Figure17. The unit vertical load is evenly imposed on the upper surface of the design field. The design domain with length and width of 3×1 is discretized by a 150×50 FE mesh. The upper bound of volume constraint is set to 0.4. A rectangle domain with length and width of 3×0.1 on the top of the design domain is a nondesign area.Figure 17
The bridge problem.The initial layout is shown in Figure18. The initial values of the design variables for each component of the QVT method and this paper method are set as 0.380.040.060.040.7 and 0.390.410.390.040.060.040.7, respectively.Figure 18
The initial designs for the bridge problem: (a) the QVT method and (b) the paper method.
(a)(b)From the intermediate steps of the optimization process shown in Figures19 and 20, the components of this paper method have more flexible deformation ability in the optimization process, and the smooth optimization structural boundaries with circular arc degree are obtained. Although the number of components in the design domain is small, the effective load transmission paths can still be obtained by proper bending deformation of the components. The objective function values are Obj=155249.20 and Obj=146920.09, respectively. It can be observed that the curved bridge structures have a smaller objective function value.Figure 19
Some intermediate optimization results and the final results of the bridge example with the QVT method.
(a)(b)(c)(d)Figure 20
Some intermediate optimization results and the final results of the bridge example with the paper method.
(a)(b)(c)(d)


## 6. Conclusions

In the present paper, a topology optimization method based on the MMC to solve the local nonsmooth problem of components connection areas is proposed. It enhances the variability of components’ shape by coordinating multiple length and thickness variables, and the smooth structural boundaries can be formed in the connection areas between the components. This method preserves the advantages of the original MMC approach and enhances its geometry modeling capabilities. The effect of undesirable positioning or deformation of some components in a layout with more components on the rate of convergence can be optimized by the hierarchical variables calling strategy. Several numerical examples are used to verify the effectiveness of the proposed method. It is also observed that the initial layout of the components is very important for the rate of convergence. In principle, the hierarchical variables calling strategy can be applied to the structural topology optimization that needs to optimize the initial layout of components to improve the rate of convergence. But how to achieve a reasonable initial components layout way is a more challenging task, and the related research will be reported in separate works.



---

*Source: 4749698-2020-10-22.xml*