Three-dimensional stress-based topology optimization using SIMP method

. Structural topology optimization problems have been formulated and solved to minimize either compliance or weight of a design domain under volume or stress constraints. The introduction of three-dimensional analysis is a more realistic approach to many applications in industry and research, but most of the developments in stress-based topology optimization are two-dimensional. This article presents an extension of two-dimensional stress-based topology optimization into three-dimensional using SIMP method. The article includes a mathematical model for three-dimensional stress-based topology optimization problems and sensitivity analysis. The article also includes ﬁ nite element analysis used to compute stress induced in the design domains. The developed model is validated using benchmark problems and the results are compared with three-dimensional compliance-based formulation. From the results, it was clear that the developed model can generate optimal topologies that can sustain applied loads under the boundary conditions de ﬁ ned.


Introduction
Structural optimization problems that deal with an assemblage of materials to sustain loads in the best way have been solved using either size, shape, or topology optimization. In the case of size and shape optimization, some predefined size and shape of solid elements are included in addition to the loading and boundary conditions before the optimization starts. Structural topology optimization is a method used for an optimal material distribution in the design domain without prior knowledge on the shape and size of the design domain. This freedom gives the designer and end users to have structures with lightweight, high performance, and aesthetically fit and safe design.
Solid Isotropic Material with Penalization (SIMP) [1,2], Homogenization Method [3,4], Level Set Method (LSM) [5,6], Evolutionary Structural Optimization (ESO) [7,8], and Bidirectional Evolutionary Structural Optimization (BESO) [9,10] are among the methods that have been used to formulate a structural topology optimization problem. Among the methods used for formulation of optimization problems, SIMP is the most common due to its simplicity and high computational efficiency [11]. Different algorithms have been suggested and implemented to solve an optimization problem, including optimality criteria method, convex linearization method, successive linearization method, and sequential linear programming. Optimality criteria method is the most common algorithm used to solve optimization algorithms [12][13][14][15].
Most of the developments and researches on structural topology optimization are formulated and solved for minimizing compliance [15][16][17]. The difficulty in considering multiple loading, absence of consideration for self-weight of a structure, and stress and displacement constraints make difficult optimal layouts from this formulation to be utilized in real-world scenario [11]. To address these limitations, different methods have been proposed to include stress constraints in the design domain for various objective functions [16,18,19].
The developments in structural optimization that uses SIMP method focus on two-dimensional problems and those that are formulated for three-dimensional problems consider a compliance minimization problem [15,20]. So far stress three-dimensional structural topology optimization problems are formulated and solved using LSM [21,22] and BESO [23]. This article aims to provide stressbased topology optimization mathematical model for three-dimensional optimization problem using SIMP method.
A stress-based topology optimization (STO) is stated to minimize the weight of a structure using SIMP method. The problem stated in equation (1) is formulated using a von Mises stress failure theory: Subjected to : g x e ð Þ ¼ s vm s yield < 1 : Here, V = the volume of structure to be minimized P e = design variable v = elemental volume P = penalization factor s vm = von Mises stress induced in each element s Y = yield strength of a material (yield stress) K = global stiffness matrix F = global force vector According to this theory, a material will fail when the von Mises stress induced in the material exceeds the yield strength of a material as shown in equation (2).
where s 11 , s 22 , and s 33 are the principal stresses along X, Y, and Z axes, respectively; s 12 , s 23 , and s 13 are the shear stresses in XY, YZ, and XZ planes, respectively.
The stress constraint is relaxed by an approach suggested by Duysinx [24] to avoid the singularity phenomenon associated with the discontinuity of the constraint function due to the removal of elements as shown in equation (3).
where j is a relaxation parameter value in the range of [0.001-0.1]. For relating the macro stress to the micro stress levels, a local stress interpolation proposed by reference [25] is used as shown in equation (4): where s(x) is the local stress of material point. The exponent q > 1 is used for preserving physical consistency in the modeling of a porous SIMP material [16]. Macroscopic elasticity tensor can be related to the constitutive elasticity tensor using power law approach as shown in equation (5).
e(x) macroscopic average strain of a material point D e (x) macroscopic elasticity tensor Here, D 0 = constitutive elasticity tensor (6): The design domain is assumed to be rectangular and discretized by brick finite elements as shown in Figure 1. The average strain of an element at the centroid of an element can be expressed as equation (7): B e is the strain displacement matrix, which will be used to calculate stiffness matrix k i (r e ) defined in equation (8), can be calculated as shown in equation (10) and u e is an elemental displacement vector. where j, h, and g are the natural coordinates along X, Y, and Z axes, as shown in Figure 1, respectively. ::: Here, N1-N8 are shape functions in natural coordinate systems defined by equation (11): The degree of freedom for an element as shown in equation (13) is used to extract elemental displacement vectors, using node numbers defined in terms of number of elements defined in X,Y, and Z, as shown in equation (12): 3n4 À 2 3n4 À 1 3n4 where nelx, nely, and nelz are number of elements in X, Y, and Z, respectively, and n1, n2, n3, and n4 are node numbers of a brick element.
The global displacement matrix where elemental displacements are extracted using an elemental degree of freedom is a solution of the equilibrium equation (14): where K is the global stiffness matrix and U and F are the global displacement and force vectors, respectively. The elemental displacements can be extracted from the global using the degree of freedom as shown in equation (15): Having the definition for strain and constitutive matrix, stress state defined in equation (4) can be expressed as equation (16): From the above derivations and relations, the stressbased topology optimization problem in equation (1) can be expressed as shown in equation (17).
The design variable is relaxed from the lower boundary to avoid the discontinuity of the stress constraints and stiffness matrix.

Subjected to:
Px P Àq e D 0 Bu e s yield < j 1 À x e ð Þ : An optimality criteria method is used to solve the defined problem, by which updating scheme for design variables can be expressed as equation (18): where m is a positive limit, which usually takes a value of 0.2,h is a numerical damping coefficient with a value of 0.5 [26,27],b is dependent on the type of problems defined as shown in equation (19): where ∂v ∂x and ∂g ∂x are derivative of objective and constraint functions with respect to the design variable, respectively, and l is a Lagrangian multiplier, where the value can be found using a bisectioning algorithm. The sensitivity analysis can be found as shown in equation (20): The derivative of the objective V(r e ) with respect to the design variable r e can be calculated as The derivative of the constraint function g(r e ) with respect to the design variable r e can be calculated as The derivative of displacement vector can be calculated from the equilibrium equation defined in equation (14) [15] and STL for optimal layout from (b) developed model and (d) [15].
Derivative of the stiffness matrix can be calculated from equation (8) as Substituting the derivative of the stiffness matrix in equation (25), the derivative of the displacement vector becomes ∂UðrÞ ∂r ¼ Àðx P k o e Þ À1 px P À1 k o e UðrÞ ð 26Þ The derivative of the constraint function can be expressed as

Numerical results
Validating the developed model benchmark problems includes simply supported beam, classical cantilever beam, and L-shape beam, which are simulated using a MATLAB code described by the flowchart in Figure 2. All the design domains are discretized by a brick finite element. The material considered for all benchmark problems described in this article has a Young's modulus of E = 1.0 MPa, a Poison's ratio of v = 0.3, a von Mises stress of 10 Mpa subjected to a unit load, and a penalization factor of 3. A three-dimensional compliance-based structural topology optimization developed by Liu and Tovar [15] is used for validation of the proposed mathematical model. The maximum stress induced in the design domains and compliance is calculated and compared during each iteration for both formulation techniques.
The design domain and loading conditions for each case study are defined followed by STL (stereo lithography) file for optimal material distribution. The maximum stress induced in both compliance-based formulation is calculated during each iteration. Compliance of a design domain for stress-based formulations is calculated during each iteration and compared with corresponding compliance-based formulations.
STL files for each optimal plot are generated using a MATLAB code [15] by using 0.5 as a default cutoff value. The developed model is compared with a three-dimensional compliance-based model developed by Liu and Tovar [15]. All benchmark problems considered in this article are simulated using MATLAB on Dell Ài5-4 GB RAM (3.41)-3.2 GHz computer.

Classical cantilever beam
The first case corresponds to a cantilever beam under loading and boundary conditions as shown in Figure 3. The design domain is discretized into 60 Â 40 Â 4 [15] in X, Y, and Z axes, respectively. Figure 4a and Figure 4c show optimal material distribution from the developed model and a compliancebased topology optimization [15], respectively, for Figure 3a. From the figures, the optimal plots, using stress-based topology optimization yields a smaller number of transition elements. The disconnected elements as shown from the equivalent STL files in Figure 4b and Figure 4d will make the manufacturing of the final optimal layout difficult. Figure 5a shows a variation of maximum stress induced in the design domain. As it can be seen from the figure, the maximum stress induced in the design domain while formulated using STO is less than the classical compliance minimization formulation. This shows that for those designs where stress is the main design variable, it is advisable to consider optimal layouts from STO-based formulations. Figure 5b shows a variation of compliance for the design domain under the two formulation techniques. It is noticeable from the figure that the difference in compliance of the design domain is less as the optimization converges.
The second case study considered is a cantilever beam with a predefined shape in the design domain as shown in Figure 3b. Figure 6a and Figure 6b show optimal material distribution from the developed model and a compliancebased topology optimization [15], respectively. The equivalent STL files are shown in Figure 6c and 6d under STO and compliance minimization formulation, respectively. The number of disconnected is more in optimal layouts from compliance minimization formulation, which increases the existence of unwanted material in the final output. Figure 7a shows variation of maximum stress induced in the design domain. As it can be seen from the figure, the difference in maximum stress induced in the design domain is less when the optimization is converged. Figure 7b shows variation of compliance for the design domain under the two formulation techniques. It is noticeable from the figure that the difference in compliance of the design domain is less as the optimization converges and the formulation based on STO takes more iteration number to coverage.

Simply supported beam
The third case study considered is a simply supported beam under the loading and boundary conditions, as shown in  [15] and STL for optimal layout from (c) developed model and (d) [15]. Figure 8. The design domain is discretized into 60 Â 20 Â 4. Figure 9a and Figure 9b show optimal material distribution from the developed model and a compliance-based topology optimization [15], respectively. As it can be seen from the figures, the optimal plots using stress-based topology optimization yields a smaller number of transition elements. The percentage of solid elements, transition, and void elements in the respective formulation techniques is 59.58, 40.34, 0.08 and 59.58, 40.25, 0.16, respectively. From the equivalent STL files in Figure 9c and Figure 9d, the optimal layout based on STO is less complex when the problem is formulated based on classical compliance minimization. Figure 10a shows variation of maximum stress induced in the design domain. From the figure, maximum stress induced is less when the design domain is formulated using STO at the cost of computational time. This shows that for those designs where stress is the main design variable, it is   [15] and STL for optimal layout from (c) developed model and (d) [15]. Fig. 10. Variation of (a) maximum stress and (b) [15] and STL for optimal layout from (c) developed model and (d) [15]. advisable to consider optimal layouts from STO-based formulations. Figure 10b shows variation of compliance for the design domain under the two formulation techniques. It is noticeable from the figure that the difference in compliance of the design domain is less as the optimization converges.
From the simulation results, it can be shown that the developed model can generate optimal material distribution for three-dimensional stress-based structural topology optimization.

Conclusion
Structural topology optimization that seeks optimal material layout has been formulated and solved based on compliance minimization. Some efforts have been done to formulate and solve, including stress constraints. Though considering stress in the formulation and solution of optimization problem is more acceptable from an engineering point of view, most of the developments are limited to two-dimensional cases.
This article presents a mathematical model for threedimensional stress-based topology optimization based on SIMP method. The developed model is tested with benchmark problems. A brick element is used for discretizing design domains and finite element analysis where stiffness matrix and displacements are calculated and used for computation of objective and constraint functions. The article also includes representation of nodal values and formulation of degrees of freedom for each element in the design domain. An optimality criteria method is used to solve all the case studies considered in this article. For validation of the developed model, the simulation results are compared with compliance minimization results. From the results, it was clear that the developed model can generate optimal topologies that can sustain applied loads under the boundary conditions defined.