Influence of Density-Based Topology Optimization Parameters on the Design of Periodic Cellular Materials

The density based topology optimization procedure represented by the SIMP (Solid isotropic material with penalization) method is the most common technique to solve material distribution optimization problems. It depends on several parameters for the solution, which in general are defined arbitrarily or based on the literature. In this work the influence of the optimization parameters applied to the design of periodic cellular materials were studied. Different filtering schemes, penalization factors, initial guesses, mesh sizes, and optimization solvers were tested. In the obtained results, it was observed that using the Method of Moving Asymptotes (MMA) can be achieved feasible convergent solutions for a large amount of parameters combinations, in comparison, to the global convergent method of moving asymptotes (GCMMA) and optimality criteria. The cases of studies showed that the most robust filtering schemes were the sensitivity average and Helmholtz partial differential equation based filter, compared to the Heaviside projection. The choice of the initial guess demonstrated to be a determining factor in the final topologies obtained.


Introduction
Periodic cellular materials (PCM) can be found in nature structures like bones, honeycombs, foams, they are also artificially produced by means of different methods for several applications such as: sandwich plates cores [1] and foams for packaging materials [2,3]. Using PCM, lighter structures are possible to design [4], and new materials with tuned mechanical properties can be created such as bulk or shear modulus maximized, negative Poisson ration and materials with an specified stiffness tensor [5,6].
Cellular materials are increasingly used in many industrial applications. In order to predict the structural behavior in the macroscopic scale, it is necessary to consider a full-scale model taking into account the internal architecture of the cellular material. If the cellular material architecture is implemented into the whole model geometry, then the computational cost for the model become costly. On the other hand, the cellular material structure effect can be taken into account in the material constitutive model using the properties obtained from the homogenization procedure of a single unit cell. To obtain the properties of PCM, it is necessary to implement some mathematical approach, in order to compute the overall material properties, these techniques are called homogenization procedures. Sigmund [5] used the homogenization method to obtain the equivalent mechanical properties for a unit cell, which represent a PCM, this approach defines the material properties in a coarse or macroscopic scale from the properties of its mesostructure or unit cell. The homogenization procedure proposed by Hassani and Hinton [7], considers the equations for obtaining the stiffness tensor by means of an asymptotic expansion, these homogenization equations are expressed as a function of the strain fields and the base material properties.
The topology optimization method was first used for periodic materials design by Sigmund [5]. The first mesostructure obtained by means of this method for continuum-like materials were based on thickness variation of the elements, instead of "artificial density" as in the latest topology optimization procedures [8].
Another popular topology optimization method is the Bidirectional Evolutionary Structural Optimization (BESO). Huang et al. [9] has implemented successfully the method BESO to obtain mesostructures that maximizes bulk and shear modulus. Additionally, the Level Set Method (LSM) was successfully used by Wang et al. [10] to propose a new design procedure for multi-phase materials with prescribed mechanical properties. Although the LSM method does not require the implementation of filters to solve the checkerboard problem, it is common to find some numerical difficulties related to the use of finite differences schemes to solve the evolution problem, for some advanced problems could be complicated to satisfy the Courant-Friedrichs-Lewi condition (CFL). Some authors have implemented the LSM method to solve the minimum compliant design problem (MCD) [11], and practical design applications such as optical and carpet cloaks [12,13].
The most popular technique for structural topology optimization is the density based approach viz. "Solid Isotropic Material with Penalization" (SIMP) developed by Bendsøe, M. P., Sigmund [8]. SIMP method is widely used to solve optimization problems, Zuo and Saitou [14] used this method to solve multi-material MCD problems and Bai and Zuo [15] developed a formulation to analyze hollow structures. Regarding PCM Xia and Breitkopf [16] proposed a Matlab code to obtain material mesostructures maximizing shear and bulk modulus. In the same work, an objective function to generate mesostructures with negative Poisson ratio was also implemented. The code developed by Xia and Breitkopf [16] was based on a Matlab code created by Andreassen et al. [17] for the solution of the minimum compliance design problem. Some other authors have used the same method to design architected cellular materials [6,18,19]. A complete review addressing the topic of topology optimization for architected materials can be found in [20,21].
Ingrassia et al. [22], Chiu et al. [23] studied the influence of optimization parameters in BESO for MCD problems. In the same way, Zuo et al. [24] assessed the performance of MMA (Method of Moving Asymptotes) and GCMMA (Global Convergent Method of Moving Asymptotes) optimization solvers to deal with MCD problems using the SIMP method. In the same work Zuo et al. [24] proposed a hybrid formulation of MMA and GCMMA to improve the efficiency and convergence of the solution. Additionally, several benchmark problems were solved using the SIMP method in references [25,26].
In the present work, we show the dependence of the solutions obtained using the SIMP method on some parameters such as: penalization factor, filtering method to obtain "0-1" (black and white solutions), optimization solvers, mesh size and initial conditions. Additionally, some issues related to the uniqueness and "monotonicity" of the objective function used in the inverse homogenization procedure for PCM were analyzed. Finally, the characteristics of the objective function and the influence of the aforementioned parameters , for inverse homogenization in plane stress problems were studied in this article.

Homogenization of Periodic Cellular Materials
To obtain the macroscopic mechanical properties of PCM, it is necessary to implement a homogenization procedure in order to estimate the effective stiffness or elasticity tensor in each iteration during the optimization process. If a representative unit cell with periodic boundary conditions is considered, the homogenized elasticity tensor in the domain Ω which represents the material macroscopic properties can be calculated using the asymptotic expansion procedure as follow [9]: where, |Y| denotes the cell volume in the domain Ω.ε kl pq are the linearly independent unit strains: ε 11 pq = (1 0 0) T ,ε 22 pq = (0 1 0) T andε 12 pq = (0 0 1) T ,ε kl pq is the induced strain field and C ijpq is the elasticity tensor for the base material. If an energy based approach is implemented, then, the unit test strains ε A(kl) pq should be induced on the boundaries which corresponds to the strain fields ε kl pq −ε kl pq . Equation (1) can be re-written in terms of element mutual energies as follow [5,16]: Equation (2) can be written on its finite element form, considering a unit cell discretized into N finite elements and the displacement solutions u A(kl) e corresponding to the strain fieldsε kl pq as follow: where, k e is the element stiffness matrix.

Elastic Problem Solution
The strain fieldε kl pq is computed from a set of three load cases (for plane stress), as showed in Figure 1, and the equations can be written as follow [18]: where, u ab is the displacements vector, e a is the unit vector along the a − th direction, 1 2 e a e T b + e b e T a are the tractions on the boundary ∂Ω, and n is a unit vector normal to the boundary.  The plane stress problem described in Equation (4) is solved using periodic boundary conditions. If it is considered a periodic displacement field, then, the induced strainε ij pq is expressed as the sum of a macroscopic displacement field and a periodic fieldũ p [16], as follow: where, y q is the length of the cell in the direction q. Considering the cell depicted in Figure 2, the displacements on the opposite boundaries are defined by: where, superscripts "k+" and "k−" define two pair of opposite parallel boundary surfaces. Eliminating u p from Equation (6) the following periodic boundary condition was obtained: this condition relates the displacements for opposite faces of the cell, giving the periodicity condition necessary to model the PCM response. The stiffness tensor is obtained by means of solving Equation (4) with the boundary conditions in Equation (7) and using the homogenization procedure in Equation (1).

Optimization Model
The SIMP method for topology optimization was implemented. The domain Ω was discretized into N finite elements and the design variables were represented by the "artificial density" of each element ρ e . Using SIMP, the Young modulus for an element E e is defined as a polynomial function [8,16]: where, E 0 is the Young's modulus for the base material, E min is a factor to avoid the singularity of the element's stiffness matrix, and p is a penalization factor, ρ e is a continuous scalar value defined in the interval 0 ≤ ρ e ≤ 1. The inverse homogenization procedure to obtain PCM can be stated as the following optimization problem [16]: where, ν e is the element volume, ϑ is the prescribed volume fraction for the PCM, K is the global stiffness matrix U A(kl) and F (kl) are the global displacements and forces vector for each load case defined in Figure 1. The authors Xia and Breitkopf [16], Chen et al. [27] and Neves et al. [6] have proposed different objective functions for Equation (9). For the maximization of shear modulus c = −(C 1212 ), bulk modulus c = −0.25(C 1111 + C 1122 + C 2211 + C 2222 ) or material with negative where β is a fixed parameter defined in the interval 0 < β < 1 and l is the current iteration during the optimization process. In the general case of Equation (9), the inverse homogenization problem formulated to obtain materials with a prescribed stiffness tensor can be stated as [18]: where, C 0 ijkl is the prescribed stiffness tensor and C H ijkl is the homogenized tensor obtained in the current iteration.

Sensitivity Analysis
In order to asses different approaches to solve the optimization problem (9), three gradient-based methods were implemented. The methods considered in this work were: Optimality Criteria [8], Method of Moving Asymptotes [28] and the Global Convergent Method of Moving Asymptotes [29]. For implementing the approaches, it is necessary to calculate the derivatives of the objective function and constraint equations with respect to ρ. The derivative of the stiffness tensor ∂C ijkl (ρ)/∂ρ was computed using the adjoint method [30] and is given by: where, k 0 is the element stiffness matrix for the solid material. The first constraint in Equation (9) is the equilibrium condition and is forced during the finite element solution. If a uniform mesh is implemented with an element volume v e = 1, thus, the volume derivative ∂V/∂ρ e = 1. The derivative of the function ||C 0 ijkl − C H ijkl (ρ)|| 2 can be obtained considering the norm definition of a fourth order tensor as the double inner product ||T|| 2 = T : T, as follow:

Filtering
Three different filtering strategies to obtain manufacturable patterns were compared. The most common filter approaches were taken into account, viz., Field averaging, Heaviside projection and filter based on the solution of the Helmholtz partial differential equation [31][32][33].

Field Averaging
In the field averaging approach the sensitivities are modified as follow: where, r min is the filter radius, N e is the set of elements i for which the distance ∆(e, i) ≤ r min and H ei = max(0, r min − ∆(e, i)) is a weight factor. In addition, the field averaging method can also be applied to ρ e , giving:ρ

Heaviside Projection Filter
A modification of Equation (14) was proposed by Guest et al. [34] . In this filter the densityρ e is projected to a physical densityρ e , where,ρ e = 1 ifρ e > 0 andρ e = 0 ifρ e = 0. The following function is introduced in order to use a gradient-based optimization scheme: where, β is a smoothing parameter. To obtain the sensitivity of Equation (15), the derivative ofρ e with respect toρ e , is computed as:

Filter Based on Helmholtz Partial Differential Equation
Lazarov and Sigmund [33] used the Helmholtz PDE solution with homogeneous Neumann boundary conditions to propose a new density filter, as follow: where, ψ is the unfiltered design field andψ is the filtered field. If a sensitivity filter is applied, then ψ = ρ ∂c ∂ρ andψ = ρ∂ c ∂ρ .

Optimization Solvers
The most common approaches for variable updating in gradient-based topology optimization procedures are: Optimality criteria (OC), Method of moving asymptotes (MMA) and Global convergent method of moving asymptotes (GCMMA) [28,29,35]. A comparison of these methods is given in [26].

Optimality Criteria
In the optimality criteria it is implemented an heuristic updating scheme defined as follow: where, m is a positive move limit, η is a damping coefficient, and B e is related to the optimality condition as function of Lagrange multiplier λ:

MMA and GCMMA Methods
MMA is based on a special type of convex separable approximation, using a first order Taylor series expansion of the objective and constraint equations. In this method, a subset of optimization problems was solved using asymptotes L j and U j as lower and upper bounds in each sub-problem, the values of the asymptotes are normally changed between the iterations. The approximation functions for each sub-problem i are chosen as: The GCMMA is an extension of MMA in which p ij are not zero simultaneously, and they are function of a non-monotonous parameter. GCMMA includes inner iterations, and the approximations used are more conservative. However, the procedure is slower than MMA. More details of these methods can be found in references [28,29].

Optimization Parameters
The solution of a topology optimization problem using the SIMP method, depends on variables such as the penalization factor p, mesh size, filtering scheme, filter radius r min , optimization solver, and the initial guess. The value of p defines the polynomial order to approximate the objective function and its derivative. For the minimum compliance design problem, the recommended value of p is considered as p ≥ 3. Good results have been reported in the literature using the minimum value p = 3 [8]. For the inverse homogenization problem p is in general defined by trial and error or using the continuation method [32] until to get a convergent 0-1 solution. p = 5 was used by Xia and Breitkopf [16] obtaining good results. Filtering schemes provides mesh independence and manufacturable designs [32]. The effects of different filters were evaluated for minimum compliance designs by Andreassen et al. [17] where the savings in memory usage and computational cost for the Helmholtz PDE based filter are highlighted. The Heaviside projection filter is used to achieve a minimum length scale in the optimized design and to promote 0-1 solutions as proposed in [34]. The effects of parameter r min is well known, and in general a small value of r min will produce topologies with "fine details" and larger values yield to coarser topologies. The solver MMA usually reaches a convergent solution in fewer iterations than optimality criteria, but 0-1 solutions are not guaranteed for inverse homogenization problems, using arbitrary initial guesses. GCMMA is an improvement of MMA introducing a non-monotonous approximation, but still requires a careful selection of the initial guess to obtain an stable, convergent and 0-1 solution. In the general inverse homogenization problem, the objective and constraint functions usually consist of non-convex and implicit form of design variables, therefore it is important to select carefully an appropriate solver for a specific problem [24].

Analytical Bounds
The analytical bounds of the objective function for bulk and shear modulus maximization, were computed to obtain reference values to compare the solutions achieved using the different combinations of optimization parameters. The theoretical achievable bounds for shear and bulk modulus maximization considering quasi-isotropic and quasi-homogeneous multi-phase materials were derived by Hashin and Shtrikman [36] and Andreassen et al. [37]: (1 − ν) , where, K H and G H are the bounds for maximum bulk and shear modulus respectively, K and G are the base material bulk and shear modulus respectively, and ν is the volume fraction.

Characterization of the Objective Function
Some characteristics of the optimization problem such as: continuity, existence, and uniqueness, have to be considered in order to study the behavior of the optimization parameters. The continuity requirement can be verified looking for the design variables values that make the stiffness matrix K of the Equation (9) not singular. the factor E min in Equation (8) is added in order to avoid the singularity of K. As for the solution uniqueness, it is evident that, the solution is not unique, because there are different unit cells that generate the same PCM (see Figure 3). The appropriate selection of a solver or updating scheme for the design variables, it depends on the characteristics of the objective function. The most relevant is the monotonicity, for instance, a gradient based-method is suitable when the objective function is convex and monotonic. There is not a formal definition of monotonicity for multi-variable functions, but the evaluation of the function behavior modifying the design variables (one at a time) it is accepted as a "monotonicity" test for topology optimization problems [38,39].  A local "monotonicity" test was performed for the bulk modulus maximization objective function, verifying 100 random elements: 40 at the top, bottom, left and right boundaries and 60 at the center of the domain. ρ e was varied in the interval 0.0 ≤ ρ e ≤ 1.0, with increments of 0.05, and considering the initial guesses (a) and (b) displayed in Figure 4. In Figure 5 are presented the results. It was observed a monotonic behavior of the objective function when an uniform initial guess (Figure 4a) was considered. If a fully random initial guess (Figure 4b) is used, thus a non-monotonic behavior is obtained for some elements located in the right boundary. The same test was performed for the shear modulus maximization objective function, but it was not obtained any element with non-monotonic behavior.    When the optimization function is "non-monotonic", thus one or more local optima points exist. To verify the local optima supposition, the objective functions for shear and bulk modulus maximization were evaluated using 100 different random initial guesses similar to Figure 4b, and the obtained topologies were classified into groups by means of the Pearson correlation coefficient ξ c [40], where, all the solutions with ξ c ≥ 0.9 were considered to be the same topology. All the optimizations were performed using the MMA solver, a mesh size of 100 × 100 elements, volume fraction ν = 0.5, penalization factor p = 5.0 and a filter radius r min = 5.0 using the sensitivity average filter. The Pearson correlation coefficient is defined as: (22) where, A and B are the solutions to be compared, and the subscripts m and n refers to the element index,Ā andB are the mean values of A and B respectively.
The results are presented in Figures 6 and 7, where the existence of multiple local optima points is illustrated in a frequency diagram. The representative topologies for some groups and its corresponding mean objective function value are displayed in the same figures. In Figures 6 and 7 it can be observed that the solutions with the best mean optimum value have a higher probability to be achieved for the shear and bulk modulus maximization, when the optimization was started from a random initial guess. Additionally it was observed that the same periodic material can be defined by different unit cells as defined in Figure 3, confirming the non-uniqueness of the solution. The objective functions for bulk and shear modulus maximization described in Section 3 are particular cases of the objective function for the general inverse homogenization problem stated in Equation (10), which is assumed more complex to solve than its particular cases. In order to compare the objective function for bulk modulus maximization and the function for general inverse homogenization, it was performed the same analysis presented in Figures 6 and 7, using the stiffness tensor C 0 1111 = C 0 2222 = 0.3232, C 0 1122 = C 0 2211 = 0.0470, C 0 1212 = 0.0276, obtained with the bulk modulus maximization objective function and a fully random initial guess. In this analysis, not all the solutions were feasible, only 61/100 runs were 0-1 topologies, according to the following criteria: 23) where, N is the total number of elements, N tol ρ is the number or elements with a density value within a tolerance tol and f is the required fraction. The results are presented in Figure 8, only feasible solutions were considered, in this test, the global optima it was not the most frequent solution, instead, it was found an scaled version of itself as the most common solution.  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29    Additionally to the initial guess, r min has an important effect in the solution (magnitude and topology), therefore, an additional study was carried out using five different filter radius in conjunction with five volume ratios, and considering the bulk modulus maximization objective function, for the initial guesses shown in Figure 4a,b, the volume fraction and filter radius were varied within the intervals (0.3 ≤ v ≤ 0.7) and (0.3 ≤ r min ≤ 5.0) respectively. The other parameters were defined using the same values implemented to obtain the topologies in Figure 7. The results are presented in Figure 9 for the initial guess in Figure 4a and in Figure 10 for the initial guess in Figure 4b. r min values that did not appear in the Figure 10 are non-convergent solutions. Additionally the Hashin and Shtrikman [36] analytic bounds are presented in Figures 9 and 10. It is observed in Figure 9 that in general the highest bulk modulus was obtained with the minimum filter radius when the iteration process started from a initial guess near to the optimum such as in Figure 4a, this is not valid for the fully random initial guess in Figure 4b, where, the best optima is not necessarily obtained with the lowest r min value as observed in Figure 10.

Interaction among Parameters
The initial guess and filter radius are not the only factors that affect the solution, in addition, the penalization power p, filter type, mesh size and updating scheme have also an important role in the solution. A fast and robust algorithm to solve topology optimization problems should be able to reach a feasible solution in a wide range of values for all the optimization parameters. The solution is considered feasible if: • It is not a Grey solution (Figure 11a,b show examples of gray solutions, which are not desirable) • The convergence is reached within a reasonable number of iterations (500 it was considered for all test) • The final topology is not disconnected (see Figure 11d) • Early convergence is not presented (see Figure 11c To analyze the interaction of all aforementioned optimization parameters, the runs analogous to a full factorial design of experiments [41] were completed, considering the factors and levels presented in the Table 1. Taking into account all the levels combinations a total of 1478 runs are obtained. Once completed all the runs, were found only 980/1478 feasible solutions, according to the feasibility criteria defined above. The obtained results are summarized in the Table 2, with the best optimal solution ID (as defined in Figure 12) for each level, and the number of feasible solutions.The most robust solver found in the analysis was the MMA, with this method, were obtained a large mount of feasible solutions when it was used in combination with all the other parameters, the best optimum value obtained with this method was achieved using the heaviside filter, a filtering radius r min = 3.0, penalization factor p = 5, a mesh of 150 × 150, an uniform grid initial guess, and the filter applied to the sensitivity field.   In order to obtain a detailed view of each solver's performance, the same results presented in Table 2 are shown in Tables 3-5 but displaying the results of each solver in combination with all the other levels.
In the Table 3 can be observed that all the early convergence solutions were obtained using the optimality criteria method, with the fully random initial guess, in fact, feasible solutions are only obtained using the optimality criteria with the center hole initial guess which is relatively close to the global optima. The MMA method proved to be the most robust, achieving convergent solutions starting from all the initial guesses, the method fails when the filter is applied to the density field instead of the sensitivity field. On the other hand, the GCMMA had problems to reach convergent solutions starting from an uniform grid initial guess, with small penalization factors and with the Heaviside filter scheme. An advantage of GCMMA is that convergent solutions can be obtained applying the filter to the density, as well as, the sensitivity field. A generalized difficulty, to get feasible solutions using the heaviside projection filter was observed. The codes implemented to obtain the results presented in this work, can be found in the Supplementary Materials.

Conclusions
The influence of optimization parameters in the solution of inverse homogenization problems using the SIMP method was studied, in addition, the features of the objective function were analyzed and it was demonstrated numerically that for some solvers and filter schemes the optimization problem become ill posed. The MMA was found as the most robust solver. Using the field average filter were obtained the largest amount of feasible solutions, and the heaviside projection filter presented problems to achieve convergent solutions in combination with the most of the studied parameters. The choice of the initial guess demonstrated to be a determining factor in the final topologies obtained.