Multi-material Topology Optimization Based on Multiple SIMP of Variable Density Method

 Abstract: In this paper a multiple SIMP (solid isotropic material with penalization model) of variable density method is proposed to solve the problem of muli-material topology optimization. All candidate materials including void material are arranged in descending order of elastic modulus. The material conversion scheme of multiple SIMP is based on the elastic modulus of the candidate material after interpolation. The iterative criterion of muti-material topology optimization is derived from the Kuhn-Tucker condition using the guide-weight method. Three examples show that it is effective and moderate to use the proposed method to solve the problem of muti-material topology optimization.


Introduction
As a kind of structural optimization, topology optimization is a mathematical method for optimizing the distribution of materials in a given design domain according to given load conditions, constraints and performance indicators [1][2][3].In the past decades, a variety of topology optimization methods such as level set method [4][5], phase field method [6][7] bi-directional evolutionary structural optimization (BESO) method [8][9], and variable density method [10][11] has been extensively and deeply researched.These topology optimization methods gradually expanded from single-material to multi-materials [12][13][14][15].Guo el al. [16] proposed a level set-based approach to solve the stress-related topology optimization problem.The absence of gray areas in the design domain is the most attractive feature of this approach during the optimization process.The approach is extended to stress-related topology optimization for muti-materials [17].In comparison with SIMP-based optimization approaches, the level set-based approach has lower solution efficiency.This means that it needs more iteration to get the optimal solution.
Tavakoli et al. [18] introduced a new algorithm to solve muti-materials topology optimization problem.The core idea of the algorithm is to divide the multi-material topology optimization problem into a sequence of binary phase topology optimization sub-problems.These sub-problems are solved in a sequential manner according to the traditional binary phase topology optimization problem.There are a large number of binary phase topology optimization sub-problems to be solved in the proposed algorithm.The computational cost of this algorithm increases dramatically.It usually takes over 1000 iterations to find an optimal solution.
To solve the muti-material topology optimization problem, Xia et al. [19] integrate the material conversion scheme of the BESO with the level set method.This method inherits the advantages of the BESO method and the level set method, and overcomes the shortcomings of the two methods.
Among the various popular interpolation schemes of variable density method, SIMP is a commonly used density-stiffness interpolation model.The interpolation model uses the element relatively density as the design variable, and the material properties are simulated by the exponential function of the element relatively density.The elastic modulus in SIMP for single material [20] can be formulated as ( ) where E(x) is the elastic modulus after interpolation; x is the element relatively density; p is penalization factor and E 0 is the elastic modulus of the material.
For the two-phase material problem, the SIMP model can be expressed as ( ) 12 ( 1) where E 1 and E 2 are the elastic modulus of the two materials, respectively.The three-phase material usually contains two material and void material.The SIMP model of three-phase material requires two design variables to describe.It can be further expressed as ) According to this law, the SIMP model m-phase material requires (m-1) design variables to describe.Blasques [21] applied the SIMP-like interpolation model to muti-phase material topology optimization problem of composite beam cross sections.It also requires (m-1) design variables.
Zuo et al. [22] proposed an ordered SIMP interpolation method to solve muti-phase material topology optimization problem.This method interpolates the elastic modulus for multiple materials with respect to the normalized density variables by introducing power functions with scaling and translation coefficients.The biggest feature of this method is that there is no need to introduce any new design variables.Yin et al. [23] proposed a new material interpolation model based on SIMP method, which is named the peak function model.This model can easily include multi-phase materials in the design by using a linear combination of normal distribution functions.The outstanding feature of this model is also that there is no need to add any new design variables.
Either the previous method for muti-material topology optimization has more iteration or additional design variables need to be added.In this paper, under the premise of without increasing design variables, we aim to solve the problem of muti-material topology optimization.

Multiple Variable Density Method
The structure is composed of m candidate materials μ j , with corresponding densities ρ j , and elastic modulus E j (j=1,2, … ,m).The candidate material μ m is the void material.Let all the m candidate materials be arranged in descending order of elastic modulus such that E 1 >E 2 >…>E m .The structure is divided into n elements and the initial material of all elements is candidate material μ 1 .
In the multiple SIMP of variable density method, the material properties are still simulated by the exponential function of the element relatively density.The elastic modulus in multiple SIMP of variable density method can be formulated as ( ) ( ) where ( ) Ex is the elastic modulus after interpolation by the jth material μj.j i x (i=1,2,…,n) is the relative density of the ith element whose material is the jth material μj.Ej is the elastic modulus of the jth material μj.
In the multiple SIMP of variable density method, the material of an element can change from one to another, as shown in Fig. 1.When an element changes from the jth material μ j to the (j+1)th material μ j +1, the elastic modulus after interpolation does not change.It is expressed as We can get Eq.( 6) from Eqs.( 4) and ( 5).
( ) ( ) The Eq. ( 7) is expressed in another form as The maximum value of , the element changes from the jth material μj to the (j+1)th material μj+1.At this time, the relative density of the ith element whose material is the (j+1)th material μj+1 is calculated as Figure 1 Elastic modulus in multiple SIMP of variable density method 3 Muti-material Topology Optimization using Multiple SIMP of Variable Density Method

The Mathematical Model
A matrix X is defined to store the relative density of elements.The number of matrix rows is the number of elements n, and the number of matrix columns is the number of candidate materials m.It is expressed as An element can only be assigned one of m materials.It means that only one element in each row is (0, 1] and the other elements are equal to zero.
A column matrix I is defined, whose elements are equal to 1.It is expressed as We are setting a matrix X which is equal to X times I.
It is expressed as For each i The mathematical model of muti-material topology optimization under the single working condition can be expressed as 12 min max ( , , , ) ( ) 0 where vector X is an n-dimensional vector of the design variables; f(X) is the objective function and g(X) is the constraint function.n is the number of finite element meshes and m is the number of candidate materials.x min and x max are the lower and upper limits of the design variable, respectively.

Guide-weight Method
In order to solve the multi-material topology optimization problem, we first need to construct the Lagrange equation as follows where λ is a Lagrange multiplier.
Based on the Kuhn-Tucker condition, the optimal value X* must be satisfied When the value of x is between the upper and lower limits of the design variable, we can obtain the following equation as We mainly discuss the problem of multi-material topology optimization under mass constraint in this paper.So f(X) represent muti-material structural performance and g(X) represent the total volume or mass of muti-material structure.
Using the guide-weight method [24][25][26][27][28], the four key formulas are defined as follows: where x .G i is the guide-weight of j i x and G is the total guide-weight.
Substituting equations ( 17) and ( 19) into equation ( 16), then we have 1, 2, , 1, 2, 3, , The above formula is the iterative criterion which can solve the problem of muti-material topology optimization taking advantage of the guide-weight method.When the muti-material topology optimization is performed, the iterative criterion can be written as (22) A step factor α is introduced into the iterative criterion.The purpose is to ensure convergence of optimization results.Considering the lower and upper limits of the design variable, the iterative criterion is finally indicated as Case(1):If ( ) (1 ) In this case, the element relative density will change and the material of the elements does not change.
, 1 1 1, or2, or3, or 1 In this case, the ith element changes from the jth material μj to the (j+1)th material μj+1.The relative density of the ith element whose material is the (j+1)th material ,1 ij x + is computed by Eq. (24).Because an element can only be assigned one of m materials.The relative density of the ith element in the jth , ij x is assigned a value of xmin as Eq.(25).

3.3Lagrange Multiplier λ
When the design variable X takes the optimal value X*, the Lagrange multiplier λ can be obtained by means of 1, 2, , Hence W i is the equivalent mass of muti-material structure.When performing a multi-material topology optimization, Lagrange multiplier λ can be calculated as where W 0 is the mass or volume constraint of muti-material structure.

The Minimum Compliance Problem of Muti-material Topology Optimization
In this section, we mainly discuss the minimum compliance problem of multi-material topology optimization by virtue of the above proposed method.The element relative density is selected as the design variable .The structural mass is as the constraint function.
It can be formulated as where C is the compliance of muti-material structure.M 0 is the initial mass of muti-material structure.M is the mass of muti-material structure after optimization.f 0 is a mass fraction.
How to get the explicit expression between structural compliance and the design variable is the key to solving multi-material topology optimization.The structural compliance C can be expressed as where U is the displacement vector.F is the load vector.K is the structural stiffness matrix.u i is the element displacement vector.k i is and the element stiffness matrix.
According to multiple SIMP of variable density method, we know that () where ki is the stiffness matrix of the ith element after penalization.k j i is the initial stiffness matrix of the ith element assigned to the jth material μj.The structural compliance C is rewritten based on multiple SIMP of variable density method.It can be formulated as Derivation of structural compliance C with respect to design variables is obtained as The load vector F has nothing to do with the relative density j i x , so the above equation is simplified as below.
From the equation KU=F, the following equation is obtained as From equations ( 35) and (36), we have Let's substitute equations ( 37) and (31) into equation (19), and then sort to get the following equation.
T () u κ u The total guide-weight G can be calculated as The mass of muti-material structure M is expressed as Derivation of the mass of muti-material structure M with respect to design variables j i x is obtained as [29] Let's substitute equations (30) and (41) into equation ( 17), then we have Equations ( 38) and (42) are plugged into equation ( 22), and then it can be written as 1, 2, , 1 2 The above formula is the iterative criterion which can solve the muti-material topology optimization problem of minimum compliance taking advantage of the guide-weight method.When the muti-material topology optimization is performed, the iterative criterion can be written as ( ) Taking into account the convergence of optimization and the limits of design variables, the iterative equation is finally expressed as Case(1):If ( ) (1 ) , 1 1 1, or2, or3, or 1 Lagrange multiplier λ can be computed by According to equations (44-48), the minimum compliance problem of muti-material topology optimization can be solved by virtue of the above proposed method.

Convergence Criteria
The relative error τ of the two adjacent optimization results is less than the given convergence accuracy τ max for five times in succession, that is, 5 consecutive times satisfying equation (49), it is considered that the muti-material topology optimization has converged.The convergence criteria can be formulated as 6 Numerical Experiments and Discussions

Bridge Structure under Mass Constraint
The design domain of the bridge structure [30] is 200×100 rectangular, as shown in Fig. 2. The design domain is divided into 200×100 elements during multi-material topology optimization analysis.The left lower endpoint ot the design domain is fully constrained.In the vertical direction of the design domain, the right lower endpoint is constrained.The vertical load F 1 is 20 N, located on the midpoint A of lower boundary.The vertical load F 2 is 10 N, located on the ¼ and ¾ of lower boundary, respectively.The other necessary parameters corresponding to these test cases are listed in Table 1.

6.1.1Muti-material topology optimization of bridge structure
Figure 3 shows the evolution of muti-material topology optimization for bridge structure #1 with three kinds of materials including void material.It can be seen that most of the elements are converted from material μ 1 to material μ 2 in the first iteration.As the number of optimization increases, the amount of material μ 1 is getting less and less and the amount of material μ 2 is getting more and more.In the 11th iteration, some elements are changed from material μ 2 to material μ 3 (void material).After the 16th iteration, the amount of material μ 3 is getting more and more and the main body of optimal topology is forming.Until muti-material topology optimization is completed, an optimal topology with 3 materials is obtained, as shown in Figure 3(f).The main body of optimal topology is composed of material μ 2 and material μ 3 .Material μ1 only appears in areas with high stress.

6.1.2The effect of mass fraction f on optimization results
Figure 4 shows the optimization curves of structural compliance C for bridge structure #1, #2 and #3 when mass fraction f is 0.2, 0.25, and 0.3, respectively.As the number of optimization increases, three curves of structural compliance increase rapidly to a maximum value and then fall smooth.The optimal topologies are obtained after 36, 33, and 19 iterations, respectively, as shown in Fig. 4. We find a law that the structural compliance and number of iterations when the optimal topology is obtained gradually becomes larger, as the mass fraction gets smaller.But the optimal topologies are similar, which can prove that the proposed method is effective and robust.

The effect of penalization factor p on optimization results
Figure 5 shows the optimization curves of structural compliance C for bridge structure #2, #4 and #5 when penalization factor p is 3, 3.5, and 4, respectively.The optimal topologies are obtained after 33, 18, and 15 iterations, respectively, as shown in Fig. 5.We find the second law that the number of iterations gradually decrease as the penalty factor p increases.This means that increasing the penalty factor can get a lower number of iterations.Similar optimal topologies can be obtained when the penalty factors are different, which again proves that the method proposed in this paper is effective and moderate.
Figure 5 Optimization curves of structural compliance when penalization factor p is 3, 3.5, and 4

6.1.4The effect of the number of candidate materials m on optimization results
Figure 6 shows the evolution of topology optimization for bridge structure #6 and #7 with four kinds of materials including void material.The main body of optimal topology with four kinds of materials is similar with the main body of optimal topology with three kinds of materials.The optimal topology for bridge structure #6 is composed of material μ 3 and material μ 4 .Materials μ 1 and μ 2 appears in areas with high stress.The optimal topology for bridge structure #7 also is composed of material μ 3 and material μ 4 .But material μ 2 only appears in areas with high stress and material μ 1 disappears.Because material μ 1 has a higher density, material μ 1 has been fully converted to material μ 2 under the condition of smaller mass fraction f.The other necessary parameters corresponding to these test cases are listed in Table 2.The optimal topologies are obtained after 43 and 48 iterations, respectively.We find that the main body of optimal topology under volume constraint is similar with optimal topology under mass constraint.Because the effect of density is not considered under volume constraint, the structural compliance under the condition of volume constraint is much greater than those under mass constraint.

Short MBB beam with real material
The short MBB beam structure is a rectangle with 200mm×100mm, as shown in Fig. 9.The left lower endpoint of the design domain is full constrained.In the vertical direction of the design domain, the right lower endpoint is constrained.The vertical load F 4 is 1000 N, located on the centre point D of lower boundary.The short MBB beam is divided into 200×100 finite element meshes.The topology optimization problem of short MBB beam with mass constraint f=0.3 is illustrated here.We use four real materials including steel, aluminium, magnesium and void material.The elastic modulus of steel, aluminium, magnesium and void material are E S =210GPa, E A =70GPa, E M =40GPa and E V =10 -6 E S , respectively.The density of steel, aluminium, magnesium and void material are ρ S =7890kg/m 3 , ρ A =2630kg/m 3 , ρ M =1740kg/m 3 and ρ V =10 -6 ρ S , respectively.Poisson's ratio v=0.3.
Figure 10 shows the evolution of muti-material topology optimization for MBB structure with 4 real materials including void material.The three different material combinations are: (A) Steel, aluminium and void; (B) Steel, magnesium and void; (C) Steel, aluminium, magnesium and void.Through a certain number of optimization calculations, we obtain the optimal topologies of the three material combinations.The main body of the three optimal topologies are similar.Steel is mainly distributed in high stress areas.Low stress area is occupied by aluminium and magnesium.But the details of the three optimal topologies are different.Comparison with the optimal topology of combination A, steel appears in the upper middle region of the optimal topology of combination B. In the optimal topology of combination C aluminium replaces magnesium in the upper middle region.At the same time, aluminium also appeared around the steel.When the topology optimization converges, the values of the mass fraction are different.Due to the different mass fractions, the minimum compliance problems of three material combinations cannot be reasonably compared.
The product of structural compliance C and mass fraction f when the convergence conditions are met are used as the unified performance index PI.It can be formulated as

PI C f =
(50) The dimension of the performance index PI is the same as the dimension of the objective function.The performance index PI can be called the objective function including constraint conditions.In this article it is referred to herein as structural compliance including mass fraction.The smaller the performance index PI, the smaller the objective function including constraint conditions.Therefore, it is of practical significance to use the performance index PI to evaluate the optimal topology of the three material combinations.
When the optimal topology is obtained, the parameters of the three material combinations are shown in the Table 3.After 39 iterations, the combination C obtained the optimal topology.At this time, the structural compliance is 137.80mm and the mass fraction f is 19.69%.The performance index PI is 27.13N.mm,which is the smallest of the three material combinations.So it can be considered that the optimal topology of combination C is the optimal topology of multi-material topology optimization.

Conclusions
The innovation of this paper is to transform the muti-material topology optimization problem into multiple single-material topology optimization problems.The material conversion scheme is based on the elastic modulus of the candidate material after interpolation.Compared with the single material SIMP method, the proposed method also has the advantages of fast convergence and less iteration.This paper is the first time that puts forward to solve the problem of muti-material topology optimization based on multiple SIMP of variable density method.It can provide a new thought for researchers to study muti-material topology optimization in the future.At present, the proposed method can only achieve one-way conversion of candidate materials in descending order of elastic modulus.The next step is to incorporate the idea of BESO method into my research.The purpose is to achieve bidirectional transformation of candidate materials.

Declaration
Elastic modulus in multiple SIMP of variable density method Optimization domain of bridge structure Evolution of topology optimization for bridge structure #1 Optimization curves of structural compliance when mass fraction f is 0.2, 0.25, and 0.3 Optimization curves of structural compliance when penalization factor p is 3, 3.5, and 4 Evolution of topology optimization for bridge structure #6 and #7

Figure 2 1
Figure 2 Optimization domain of bridge structureTable 1 Parameters corresponding to bridge structure Test Case Candidate materials m

Figure 3
Figure 3 Evolution of topology optimization for bridge structure #1

Figure 4
Figure 4 Optimization curves of structural compliance when mass fraction f is 0.2, 0.25, and 0.3

Figure 6 7 6. 2
Figure 6 Evolution of topology optimization for bridge structure #6 and #7 6.2 Cantilever beam structure under the condition of volume constraint or mass constraint The short cantilever beam structure is a rectangle with a length of 200 and a width of 100, as shown in Fig.7.The left boundary of design domain is fixed.A vertical load of F 3 = 10 N is applied to the centre point B on the right boundary of the design domain.The short cantilever beam structure is divided into 200×100 finite element meshes.The other necessary parameters corresponding to these test cases are listed in Table2.

Figure 7 Table 2 3 Figure 8
Figure 7 Optimization domain of cantilever beam structureTable 2 Parameters corresponding to cantilever beam structure Test Case Candidate materials m

Figure 9
Figure 9 Optimization domain of short MBB beam structure

Figure 10
Figure 10 Evolution of muti-material topology optimization for MBB beam

Table 3
Parameters when the optimal topology is obtained