Multiple materials layout optimization in a layered structure

– To obtain the optimal layout of multiple materials in a layered continuum with clear interfaces between materials, the series solid isotropic material with penalization (S-SIMP) method is proposed in this work. In practical engineering, many types of materials could be used simultaneously in a structure for reducing either cost or weight of structure. For the convenience of manufacturing, the inter-material boundaries should be simple and clear. Supposing m types of material in a structure to be layout, in the present method, in each sub-optimization, there are only k ( <m ) types of materials selected in layout optimization. Therefore, there are ( m − k +1) sub-optimizations in ﬁnding the ﬁnal layout of materials. In the ﬁrst sub-optimization, these k types of materials have higher modulus than the rest materials. Simultaneously, the rest materials are replaced with void or the softest material. After each sub-optimization, the ﬁrst material (with the highest modulus among the current materials) will be ﬁxed in the next sub-optimization. The design domain of structure decreases step by step. Numerical examples are given to demonstrate the validity and eﬃciency of the presented approach.


Introduction
With the development of high-performance computing, structural optimization methods, topology optimization is becoming more and more popular in structural design field around the world [1] over the last three decades. Several approaches were presented for finding solutions of topology optimization. Typically, material interpolation methods with homogenization method [2], variable density method [3,4], level-set method [5][6][7], evolutionary structural optimization (ESO) [8], etc. were investigated.
Research on variable density method is mature and has been applied in many structural topology optimization problems. Especially, the SIMP method [3,4] and Rational Approximation of Material Properties (RAMP)] [9] method are very popular in design engineering. In such methods, the original 0/1 design problem is relaxed into a problem with continuous design variables. But "grey material" phenomenon appears in solution. To reduce the amount of "grey" material with mid-density in materials, the artificial power law approach is suggested. The equivalent modulus in the original solid/void layout problem can be described as: where ρ ∈ [0, 1] is the design variable of density, E 1 and E 2 are the material modulus of the two phases, and typically let penalization parameter p = 3. For the "void" phase, let E 2 = 1.0 × 10 −6 E 1 [10] to avoid the singularity of global stiffness matrix in structural deformation analysis.
Similarly for the three-phase material (with two solid materials and one void) topology optimization was considered and the power-law interpolation scheme [11] is as: Where ρ 1 , ρ 2 ∈ [0, 1] are design variables in an element. For this problem, the number of design variables is doubled as comparing with the two-phase material interpolation (Eq. (1)). For example, Yin and Ananthasuresh [12] formulated a compliant mechanisms design model with three-phase (two solids and with void, E 3 ≈ 0) material by using a peak function material interpolation scheme. If materials number is greater than 3, e.g., marking "k" as the number of materials in structure, the equivalent elastic modulus of mixture of k materials can be expressed as: Article published by EDP Sciences And related multiple material layout optimization can be considered. Although multiple materials layout problem was studied for over two decades, it is still a challenging task in real practical structural design as considering the manufacturing. Besides, the complexity of interpolation model [13][14][15] exists during application, too. Saxena [16] adopted genetic algorithm in topology optimization of truss structure with multiple-material. Mei and Wang [17] presented an algorithm based on the level set method on multiphase topology optimization. Further, Wang and Wang [18] proposed a "color" level set method for multiphase material layout optimization. In their method, m level-set functions are required to represent a structure with n = 2 m different material phases. Zhou and Wang [20] suggested a generalized Cahn-Hilliard method for solving this problem. Stegmann and Lund [21] presented a discrete material optimization (DMO) approach to solve the orientation problem of orthotropic materials and the multiple material selection in structural design. Han and Lee [22] suggested a material-mixing method based on ESO for multi-phases material topology optimization. Computing 'pseudo-sensitivities' of objective and constraint functions on discrete design variables, Ramani [23] gave a new multiple material layout optimization method.
Layered structure is common engineering. Theoretically, the multiple elastic materials layout optimization can be solved by the approaches introduced above. However, the final material layout should meet the needs of practical application on consideration of manufacturing techniques. For example, the interface among many materials should be simple and clear for the convenience of manufacturing. Wang and Zhou [24] presented a phasefield method to find multiple materials layout. In 2013, Allaire et al. [19] studied the interfaces among multiple materials in an elastic continuum structure by the level set method. Considering the characteristics of the level set method, the interface between two materials is obtained by smooth interpolation. Hence, the inner boundaries of materials in their numerical examples look quite clear.
Clearly, if the number of materials is small, e.g., k = 2 or 3, the interfaces are simple and easily found. In the present study, based on SIMP method, we suggest another approach to give multiple material layout optimization with simple and inner interfaces. The method is called the series SIMP (S-SIMP) because the whole multiple material layout optimization is replaced with a series of submultiple material layout optimization (simply called as sub-optimization) based on SIMP method. In each suboptimization, there are only few materials considered as design materials. At each end of sub-optimization, the boundary of material with highest modulus is fixed and the related area will be considered as non-design domain in the subsequent sub-optimizations.
In Section 2, the description of the S-SIMP method is presented according to a stiffness design of a structure with multiple materials. The sensitivity analysis is also demonstrated. Some numerical examples are given in Section 3. Finally, some conclusions are drawn in Section 4.

Elastic design domain
In this study, the materials in structure are linear elastic and small deformation is considered. In a design domain Ω ⊆ R 2 or R 3 with Γ 1 as Dirichlet boundary condition, and Γ 2 as Neumann boundary condition (and Γ 1 ∪ Γ 2 = Γ ), the control equations are written: At where f is the body force applied on the structure, u denotes the displacement field, u 0 is the given displacement on Γ 1 , T is the boundary traction force applied on Γ 2 . σ is the stress tensor, ε is the strain tensor and D denotes the elasticity tensor.

General optimization model for multiple materials layout
Commonly, the stiffness design of a continuum with multiple materials can be expressed by the following formula when the deformation of structure is solved by finite element method.
the relative densities (design variables) of the e-th element in design domain. "m" is the total number of materials in the whole design domain. The objective is the mean structural compliance. F , U and K are the global force vector, global displacement vector and global stiffness matrix of structure. g j (j = 1, 2, . . . , J) are the related constraints considered in design process, e.g., volume constraints on the m types of material, displacement of a point on boundary, etc. J is the number of constraints in optimization.
If we only consider the optimal layout of a structure with specified volumes on materials, the constraint functions in equation (5) are written: where "EL" means the total number of elements in the design domain.

Series SIMP method
In series SIMP method, the original m types of multiple material layout optimization will be replaced with a series of k types of multiple material sub-optimization. Hence, the classification of materials for each suboptimization should be carried out before operation.

Classification of materials in sub-optimization
In stiffness design, the task is to find the optimal load transfer paths (LTP) on which the stiffer material will be layout. If there are many materials to be layout, the stiffest material should be layout on the main path to contribute more stiffness for the whole structure. Therefore, we firstly name the materials orderly from "M 1 " with the highest modulus to "M m " with the lowest modulus (see Fig. 1a). In each sub-optimization, we choose only k types of stiffest materials from the materials in design domain. Meanwhile, in the present final optimal design, the area only filled with the present stiffest material is taken as non-design domain in the next sub-optimization. Therefore, there are totally m − k + 1 sub-optimizations till finding all of the m materials layout (Fig. 1b).
In each sub-optimization, the rest materials (softer materials) will be defined as either void or the weaker material among the present selected (i.e., k types of) materials. If the summation of the volume ratios of the present selected materials is high, the "void" scheme is suggested (see Fig. 1b). On the contrary, the "solid" scheme should be adopted, in which the rest materials will be replaced with the "current" weaker material (see Fig. 1c). The "void" scheme is adopted in the present study.
Hence, if k = m, there is only one sub-optimization, i.e., it is a common multiple material layout optimization. If k = 2, it is a sequent single-material layout problem. However, if k = 2, the optimal solution of multiple material layout might be far away from that when k = m. If k > 4, the interfaces among materials will be complicated although the solution is close to the original multiple materials layout solution. Therefore, k = 2, 3 or 4 will be discussed in the following examples.

Optimization model for S-SIMP method
If k = 2, there is only one design variable in each element in design domain. Each sub-optimization can be solved by common SIMP method. If k = 3, there are two design variables, e.g., the volume ratios of the first and the second materials, the optimization model are written: where ρ i, e , ρ i+1, e are the design variables in the eth element, V 0 is the volume of design domain, f i and f i+1 are the volume fractions of the ith and (i + 1)th materials in the "current" design domain.

Sensitivity analysis
(a) If k = 2, the sensitivity of objective function is For the volume constraint The sensitivity is (b) If k = 3, the stiffness matrix of the eth element at the ith sub-optimization is Hence, ∂k So, we have When k = 3, the constraint equations are ⎧ ⎪ ⎪ ⎨ And the sensitivities of constraints are The sensitivities of objective function and constraints with respect to k > 3 can be obtained similarly.
(c) When k = 4, the stiffness matrix of the eth element at the ith sub-optimization is Hence, When k = 4, the volume constraints in the i-th suboptimization are And the sensitivities of constraints are See equation (22) next.
The sensitivities of objective function and constraints with respect to k > 4 can be obtained similarly.

Flowchart of S-SIMP algorithm
In optimization, structural deformation is numerically analyzed using finite element method. The design variables are updated by the Method of Moving Asymptotes (MMA) [25].
Step 1. Classify the materials according to their moduli, build finite element model of structure, initiate algorithm parameters, such as m, k, volume ratios of materials, and let i=1 for the first sub-optimization.
Step 4. Update the design variables in the present suboptimization by MMA method.
Step 5. Judge, if the present optimal layout of k types of materials obtained, let the elements with M i as non-design element and go to Step 6, otherwise, go to Step 3.
Step 6. Judge, if i<m-2, let i = i+1 and go to Step 2, otherwise, go to Step 7.
Step 7. Stop for post processing.
The filter scheme to overcome mesh-dependency in the present algorithm is the same as that in the work [26] by modifying the element sensitivities. The convergence of algorithm (in Step 5) is determined by two factors. One is the tolerance of design variables max ρ s − ρ s−1 η 1 (23) η 1 is a small positive integer, we take it 0.01. s and s-1 denotes the two consecutive iterations in a sub-optimization. The other is the maximum iteration number. Here we set the number to be 100. The present sub-optimization stops when any one of the two conditions reaches.

Examples and discussions
To demonstrate the validity of the present method, numerical examples with multiple materials layout are given. The objective is to minimize the compliance of the structure with volume constraints on the materials. By the way, the void scheme will be used in optimization. Therefore, the volume ratio of the stiffest material is no less than 10% in the following example.

Materials distribution in a deep cantilever beam
In Figure 2, the design domain is a deep cantilever beam with thickness of 0.01 m. The left side is fixed, on the center of the right side a concentrated force 100 N is applied. The structure is discretized into 50 × 130 elements. The Poisson's ratios are the same, e.g., 0.3. Two cases are considered.

Results of case (a)
As the optimum single-material layout is a symmetric two-bar structure, and the angle between a bar and the fixed (vertical) side of design domain is 45 • . It can be given analytically. Hence, this example is to demonstrate the validity of the present method.
In the example, we calculate the optimal layouts for k = 2, 3 and 4, respectively, to show a quantitative comparison. There are four phase materials including void. Hence, m = 4. There are three sub-optimizations for k = 2, two sub-optimizations for k = 3 and only one optimization for k = m = 4 to achieve the final optimal layout. During the process of optimization, the highest strength material is distributed at first, the rest materials layout step by step from the higher modulus till void. In Figure 3, the three final material layouts are 404-page 5   (Fig. 4). So, we find that the convergence of algorithm for k = 4 is too poor as comparing with that for k = 2 or 3. And the inter-material boundaries, especially the boundary of material (M 2 in Fig. 3), is very clear when k = 2 among the three results. Besides, the value of objective function, i.e., MSC, for k = 4 is higher than the two others for k = 2 and 3.

Results of case (b)
In Figure 5, one can find the materials layout gradient with respect to their moduli. The stiffest material (M 1) is covered by M 2. And M 3 contacts with only M 2 and M 4. Like the distributions of biological tissues in an arm of an animal, the bone which has the stiffest modulus was covered by muscle and muscle covered by skin which has the weakest modulus. Therefore, for finding the optimal layer-by-layer distributions of multiple materials in a structure, one can adopt the present S-SIMP method in design.

Medium cantilever using three materials and void
The structure shown in Figure 6 is a cantilever with sizes of 0.2 m × 0.1 m ×0.01 m. The left side of structure is fixed and a concentrated vertical load of P = 1000 N is applied on the lower right corner. The finite element mesh of structure is 100 × 50 elements.
There are two solids and void in the structure. The moduli are 2, 1 and 1.0 × 10 −6 GPa for M 1, M 2 and M 3 (void), respectively. The Poisson's ratios are the same, e.g., 0.3. The volume ratios are 0.2, 0.1 and 0.7, respectively.
In Figure 7, the two layouts of materials in structure are different from each other. Comparing with Figure 7b

L-beam using two materials and void
In Figure 9, an L-shaped structure is shown [1] with thickness of 0.01 m. The top side is fixed and a concentrated force, P = 100 N, is applied on the center of the lower right side. The structure is discretized into 25 600 elements. Three cases are considered, i.e.,  Figure 10 demonstrates the three results by different methods. Figure 10a is obtained by using S-SIMP method with k = 2. Clearly, the interfaces between materials in this figure are clearer than those in Figures 10b and 10c. Meanwhile, the layouts of M 1 in Figures 10a and 10c are very similar.

Results of case (a)
To give a further comparison between SIMP and S-SIMP methods, the modulus of material 2, E 2 , is set to be 0.2, 0.5 and 0.8 GPa, respectively. The other parameters in model keep unchanged. The optimal material layouts are shown in the Table 1.
For example, using S-SIMP method with k = 2, the layouts of M 1 in the first row in Table 1  The second row in Table 1 shows the final multiple materials layouts by using SIMP method (k = m = 3), both M 1 and M 2 are layout differently with different moduli of material 2. Without comparing their stiffness, intuitively, the material layouts are more complicate than those in the first row. Table 2 lists the MSC of structure with different moduli of M 2. When the modulus of M 2 is 0.2 GPa or lower, the structure obtained by S-SIMP method is stiffer than that by using SIMP method. However, when the modulus of M 2 is 0.5 GPa or higher, the compliance of structure designed by SIMP method is lower. This is actually caused by the difference of M 1 in structure. Figure 11 shows the final material layouts in L-beam structure with four materials (three materials plus void) by using S-SIMP method with different numbers of suboptimizations. The parameters of materials used are given in the Table 3. The final material distributions in Figure 11 are obviously different by using different methods. Figure 11a gives a layered structure. The final material distributions in Figures 11b and 11c are different slightly. After 246 iterations, the structural compliance reaches 0.402 as using S-SIMP method with k = 2; when k = 3, the compliance reaches 0.397 N.m after 182 iterations. By SIMP method (k = m = 4), the compliance approaches 0.397 N.m after 201 iterations. From the result, we can find that the compliance difference is very small between the results by S-SIMP and SIMP methods. It demonstrates that the higher moduli difference between the first material (M 1) and the rest materials leads to the optimal frame from M 1 by S-SIMP method and further leads to a stiffer structure when the amount of M 1 is not too low. If the amount of M 1 is low, the final materials distribution by S-SIMP method with k = 3 will be better than that with k = 2.

Results of case (c)
It is found that the materials in the L-beam show layer-by-layer covering the frame of M 1. The interfaces between two solids are very clear. Thus, using deposition or adhesion method, the structure can be easily manufactured.

Conclusions
In this study, a series SIMP method for multiple materials (m types) layout is proposed to find the optimal layout of materials in continuum. The original optimization is replaced with a series multiple materials layout optimization with few materials (k = types). If k = m, the optimal result by traditional SIMP method can be found. Cases E 2 = 0.2 GPa E 2 = 0.5 GPa E 2 = 0.8 GPa   However, the interfaces between materials in the structure is complicated and results in the difficulty in manufacturing. Results show that if using S-SIMP method with k = 2 or 3, the drawback can be avoided, i.e., the interface between materials is simple and clear. By the way, the computational cost by S-SIMP method is not higher than that by SIMP method. It should be addressed that using S-SIMP method with k = 2, the optimal result can be found only when the amount of the first material (which has the highest modulus) is high enough to form a simple continuous frame. Otherwise, the S-SIMP method with k = 3 should be used. Although the S-SIMP method with k = 4 can also be used to find the final materials layout, the interfaces between materials become complex and the computational cost increases.
Simultaneously, using S-SIMP method, the materials layout layer-by-layer with decreasing the modulus of material, i.e., the material i + 1 coats on the material i, will be easily manufactured in practical engineering.