A NON-GRADIENT APPROACH FOR THREE DIMENSIONAL TOPOLOGY OPTIMIZATION

. The present work is devoted to the extension of the non-gradient approach, namely proportional topology optimization (PTO), for compliance minimization of three-dimensional (3D) structures. Two schemes of material interpolation within the framework of the solid isotropic material with penalization (SIMP), i.e. the power function and the logistic function are analyzed. Through a comparative study, the efficiency of the logistic-type interpolation scheme is highlighted. Since no sensitivity is involved in the approach, a density filter is applied instead of sensitivity filter to avoid checkerboard issue.


INTRODUCTION
Topology optimization is a well-known design tool, in which the material is systematically removed from the design domain in order to get a structure with desirable properties, under specified loads and boundary conditions. For example, a minimum compliance problem searches for a design with high strengthtoweight, while a heat conduction problem is intentionally used to find structures with efficiency in conductive heat transfer. Since the introduction of 99line Matlab code by Sigmund [1], which is based on the method of solid isotropic material with penalization (SIMP), topology optimization has been drawing much attention from researchers all over the world. Andreassen et al. [2] successfully improved the efficiency of the code and reduced the total lines to 88, while keeping its simplicity and readability. The work of PolyTop was later introduced by Talischi et al. [3] for minimum compliance problems, incorporating the high accuracy of polygonal finite element analysis and the simplicity of SIMP. Other than SIMP, the level set method [4,5,6] and phase-field method [7,8,9] have also been proposed. Recently, Chen et al. [10] presented an extension of topology optimization for geometrically nonlinear structures. Regarding topology optimization involving material nonlinearity, various contributions have been made, e.g. for hyperelastic structures [11,12] and elasto-plastic behaviors [13,14,15].
Most of the works on topology optimization are based on sensitivity analysis, which involves the derivatives of objective function with respect to design variables. Such approaches are also usually referred to as gradient-based. Though sensitivity contains the physical information of the problem and thus is useful in optimization, its derivation and computation can become a burdensome task which is difficult to be implemented, e.g. in problems considering nonlinearities [16] or change of interfaces between different media [17]. Furthermore, as a mathematical aspect, the gradient-based optimization could possibly fall into local optima, e.g. optimized results may differ due to the selection of starting point [18,19,20]. In contrast, the non-gradient approaches rely only on the calculation of the objective function and do not require gradient information. The search for optimized design is usually done globally in the whole design space using some stochastic algorithm. For example, Wu and Tseng [18] employed the modified binary differential evolution, which is motivated from the evolution of species. The binary particle swarm algorithm in the work of Luh et al. [19] imitates the social behavior of animals when they collaborate to search for food. Guirguis and Aly [20] combined level set method and pattern search algorithm to develop a derivative-free tool for topology optimization. Critically, the global search of non-gradient approaches does not guarantee a global optimum, as pointed out in some works [19,21]. Another drawback is the large number of function evaluation during the stochastic search, which quickly scales up corresponding to the problem size [20,21]. Therefore, both gradient-and non-gradient approaches have their own strengths and weaknesses. In order to select between gradient-and non-gradient topology optimization, one should consider the trade-off in terms of computational complexity and efficiency.
Recently, a non-gradient technique namely Proportional Topology Optimization (PTO) has been proposed by Biyikli and To [22]. The main point of the method is that material is distributed proportionally to the compliance value calculated at each finite element. PTO was then further developed for multi-material optimization by Cui et al. [23]. However, both [22] and [23] are limited to 2D elastic problems. 3D topology optimization based on sensitivity analysis has been introduced by Liu and Tovar [24] using SIMP. Nevertheless, the number of works on non-gradient is still limited in the literature. In the present paper, PTO is extended to 3D topology optimization, for both minimum compliance and heat conduction problem types. Furthermore, a comparison is made between the power function traditionally used in SIMP for material interpolation and the logistic function proposed in [23].
The paper is organized as follows. Immediately after the Introduction, a brief review of topology optimization is presented in Section 2. Section 3 is reserved for PTO and numerical aspects to extend the algorithm from 2D to 3D domains. In Section 4, numerical examples are presented and discussed in details, demonstrating the algorithm. Comparison between the two models of material interpolation, i.e. the power function and the logistic function is also analyzed. Finally, conclusions and remarks are given in the last Section.

TOPOLOGY OPTIMIZATION FOR MINIMUM COMPLIANCE PROBLEMS
In the minimum compliance problem, the compliance (an inverse measure of structural stiffness) is minimized, while the total mass must not exceed a target value. The optimization problem reads such that: in the above equations, C is the compliance; u is the vector of displacements; K is the stiffness matrix; F is the load vector; subscript e refers to element e and ne is the number of elements. In constraint (2b), ve is the volume of element e and M is the prescribed total amount of material. It is noted that though C is in fact elastic energy, the term "compliance" has been commonly used in the research community. Linear elastic behavior and small deformation are assumed.
The design variables are elemental relative densities, e, which take value from 0 to 1, indicating whether an element is solid ( e   1 ) or voided region ( e   0 ). Applying the SIMP method, the elastic modulus is assumed to vary with respect to density as in Equation (3), E0 is the material elastic modulus and κ is a positive small number (in this work, κ = 10 -9 ) to avoid singularity when e tends to zero. Here, a power law is used to penalize the contribution of the intermediate density ( e   01 ) on the structural stiffness. Due to the penalization, the intermediate density quickly tends to 0 or 1 after some iterations. A small value of p may slow down the solution process due to insufficient penalization; while a large value of p also penalizes relatively high density (e.g. . e   05) to zero, leading to a pre-mature solution [23,24]. Typically, the penalty power, p, is chosen as p = 3. Another type of material interpolation function, recently proposed in [23], is the S-shaped logistic-type function, such that the density-stiffness interpolation reads Curious readers can refer to [23] for a parametric study of parameters b and h. In the current work, the values b = 12 and h = 3000 are adopted. It is noted that when e  = 1 , function f2 is not equal to unity. A comparison between the efficiency of the two schemes is conducted in the present work. For that purpose, here function f2 is further multiplied by a scaling factor, i.e. e h e f e he The curves of function f1 (power law) and function f2 (logistic-type law) are plotted in Figure 1 for illustration. It is observed that the contribution on the stiffness of low densities ( . 05) is small. The low-density region will be quickly penalized to 0 during the iterative loop. In power law, even the relatively high-density ( .. e   0 5 0 75 ) has a small contribution to structural stiffness and will eventually be penalized due to its insignificance, leading to the socalled "premature deletion" of the involved elements during the search for an optimized solution [23,25]. On the other hand, in the logistic-type law, a more significant number of relatively high-density are enabled to participate in the optimization process.

PROPORTIONAL TOPOLOGY OPTIMIZATION (PTO) ALGORITHM FOR 3D PROBLEMS
Originally proposed by Biyikli and To [22] for 2D problems, PTO is a non-gradient algorithm with a simple idea: the material is distributed to elements proportionally to the compliance. That is, larger values of design variables (elemental densities) will be assigned to elements that have high compliance (i.e. weak regions) and smaller values of design variables will be assigned to elements that have low compliance (i.e. strong regions). No sensitivities calculation is involved and thus the complication associated with sensitivities can be avoided. The mass constraint is directly enforced the total amount of materials that are provided, i.e. the constraint is globally imposed on the entire system. Flowchart of the topology optimization procedure is presented in Figure 2, and that of PTO algorithm is in Figure 3. Convergence of the topology optimization procedure is defined as the change between the previous vector of design variables (densities) and the current one is less than a given tolerance. After each iteration, the elemental density is updated as follows where previous e  is the density in previous iteration and optimized e  is the density obtained from PTO algorithm in the current iteration. The history coefficient α is chosen between 0 and 1 to represent the contribution of previous data to the current density. For example, α = 0 means that the density in the previous iteration does not affect the current density; and α = 0.5 means that the current density is blended equally between the previous value and the value obtained from PTO algorithm. Higher values of α means more contribution from previous solution. Thus, slower evolution can be expected. With α = 1, the previous solution is exactly kept, and therefore no evolution occurs. On the other hand, small values of α will put more weights on the optimized solution. The case α = 0 indicates that the current solution does not depend on history data at all. Stability issue (difficulty in convergence) may occur, if the optimized solution is much different from the previous solution, which likely happens especially in the first few iterations. In this paper, α = 0.5 is selected.
In PTO algorithm, TM is the total amount of material, which is determined by mass constraint in Equation (2b); RM is the remaining amount of material after distribution. If RM is still larger than a given threshold, it will then be again distributed to the elements. The process is looped until RM is smaller than the threshold, which is chosen to be 0.001 in this work. Each element receives an amount of material which is proportional to its value of compliance, i.e.
here, the elemental compliance Ce is evaluated at the centroid of each element. Notation vi denotes the volume of element i. The PTO algorithm for updating densities, as presented in Figure 3, will stop either when the remaining amount of material (RM) is smaller than a specified value or the maximum number of iteration is exceeded. If the latter is the case, the whole numerical process will be aborted due to "non-convergence". A non-convergence situation may indicate that the constraint is not appropriate, e.g. too low volume fraction. A review on the evolution of topology would be helpful to analyze why convergence does not occur.
A commonly known issue in topology with SIMP model for material interpolation is the checkerboard phenomenon, which is in fact numerical instability. To overcome this issue, some restriction must be introduced [26,27,28]. In general, a filter technique is used, either on the sensitivities or the densities [2,29]. For some certain cases, as reported in [30], the density filter may have better performance than the sensitivity one. In the current work, a non-gradient approach is employed, in which sensitivity is not calculated. Therefore, density filter is chosen in this paper. Applying filter, the density value in each element e is modified as follows ê in which Hej is the distance between the centroid of element e and element j. J is the set of elements satisfying that Hej is smaller than a given threshold Hmin. In practice, the values of Hej and the set Je can be calculated only once and stored for any element e, right after the finite element mesh is known. Extension of PTO to 3D problems brings the algorithm one step closer to practical applications (though it is currently still limited to linear elasticity and small deformation). Regarding implementation, the procedure is in general straightforward. However, it is obvious that computational cost of a 3D problem is usually much higher than a 2D problem. For every iteration, the stiffness matrix K has to be assembled and then Equation (2a) is solved for displacement. Calculation of K requires the determination of element stiffness Ke for all element e. Therefore, efficient evaluation of Ke is essential to increase overall efficiency.
It is assumed here that a finite element package for 3D linear elasticity using 8-nodehexahedral element is already available. With elastic modulus being function of density, i.e.  (4)), the constitutive matrix for an element e (isotropic material) is calculated as where D0 is the usual constitutive matrix in 3D linear elasticity in Equation (10), E0 is the material elastic modulus and ν is the Poisson's ratio. The stiffness matrix of an element e is then calculated by where Ke0 is the standard element stiffness matrix. In practice, Ke0 can be calculated only once and stored for every element.

NUMERICAL EXAMPLES
This section is served to demonstrate the extension of PTO algorithm for 3D problems of minimum compliance. For simplicity, only the eight-node hexahedral element is used in the current work. Actually, any existing finite element framework can be employed. The compliance then is calculated as a post-process of finite element analysis.
Furthermore, a comparison study is conducted for the two material interpolation schemes, i.e. the power law (Equation (3)) and the logistic-type law (Equation (5)). For that purpose, an "equivalent compliance" is defined, which is in fact the compliance value calculated using the power law for material interpolation, after the optimized densities are obtained. The optimization problem is considered to be "converged" if the change between the current solution and the previous one is less than 0.01.  Three numerical examples are considered in this section: the cantilever beam, the hollow beam, and the A-shaped bracket structure. Geometry and boundary conditions of the three problems are given in Figure 4. The cylindrical hole in the hollow beam problem is modeled by requiring elemental density values within the hole region to be zero. On the other hand, the first layer of elements on the upper surface and the first layer of elements on the clamped surface are both required to be unity, i.e. solid material.
Without loss of generality, an artificial material with elastic modulus E = 1 MPa and Poisson ratio ν = 0.3 is taken. The volume fraction, i.e. the ratio of volume of the optimized design and that of the initial model, is chosen as 0.5 for the problems of cantilever beam and hollow beam, while in the A-shaped bracket problem, the volume fraction is 0.15. Three numerical examples are considered in this section: the cantilever beam, the hollow beam, and the A-shaped bracket structure. Geometry and boundary conditions of the three problems are given in Figure 4. The cylindrical hole in the hollow beam problem is modeled by requiring elemental density values within the hole region to be zero. On the other hand, the first layer of elements on the upper surface and the first layer of elements on the clamped surface are both required to be unity, i.e. solid material.
Without loss of generality, an artificial material with elastic modulus E = 1 MPa and Poisson ratio ν = 0.3 is taken. The volume fraction, i.e. the ratio of volume of the optimized design and that of the initial model, is chosen as 0.5 for the problems of cantilever beam and hollow beam, while in the A-shaped bracket problem, the volume fraction is 0.15. Table 1, in all of three examples, the logistic-type scheme tends to need more iterations than the power law scheme to reach convergence, but lower compliance (i.e. better solution) can be obtained. This observation could be attributed to the property of the logistic-type scheme, such that it does not prematurely exclude the relatively high-density element from the optimization process. The shape of optimized structures obtained by the two schemes are presented in Figure  5, Figure 6 and Figure 7.

CONCLUSION AND OUTLOOKS
The non-gradient proportional topology Optimization (PTO) [22] has been successfully extended to three-dimensional problems of minimum compliance. Among the non-gradient techniques, the idea of the PTO is simple yet physical, i.e. material is distributed to each element proportionally to its value of compliance. As a constraint defines the total amount of material, the optimized can be obtained after some certain number of iterations. No calculation of sensitivity, i.e. derivatives of the objective function with respect to the design variables, is involved.
Comparison between the two material interpolation schemes, i.e. the well-known power law scheme and the logistic-type one, is studied in this research. It is exhibited that the logistictype scheme tends to require more iterative steps to be converged but the obtained results are better optimized than the power law scheme.
In optimization, the minimization problem and the maximization problem are equivalent. Instead of minimizing compliance, the problem can be re-formulated as maximizing stiffness. Since C in Equation (1) is a measure of compliance (in the sense that smaller C means stiffer structure), J = 1/Ccan be a representation of stiffness. Using the PTO, more material will be distributed to the elements with lower values of J.
In engineering problems, it is more interested in minimizing the structural weight while the stress value does not exceed a certain level, namely the stress-constrained problems. Sensitivity analysis for stress-constrained problems is a challenging task [31,32]. Thus, it is promising for