Issue 
Int. J. Simul. Multidisci. Des. Optim.
Volume 10, 2019



Article Number  A1  
Number of page(s)  10  
DOI  https://doi.org/10.1051/smdo/2019005  
Published online  12 March 2019 
Research Article
Threedimensional stressbased topology optimization using SIMP method
^{1}
Bahir Dar Institute of Technology, Bahir Dar, Ethiopia
^{2}
Addis Ababa Institute of Technology, Addis Ababa, Ethiopia
^{3}
Universiti Teknologi PETRONAS, Perak, Malaysia
^{*} email: hshimels278@gmail.com
Received:
24
November
2018
Accepted:
9
February
2019
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 threedimensional analysis is a more realistic approach to many applications in industry and research, but most of the developments in stressbased topology optimization are twodimensional. This article presents an extension of twodimensional stressbased topology optimization into threedimensional using SIMP method. The article includes a mathematical model for threedimensional stressbased topology optimization problems and sensitivity analysis. The article also includes finite 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 threedimensional compliancebased formulation. From the results, it was clear that the developed model can generate optimal topologies that can sustain applied loads under the boundary conditions defined.
Key words: Topology optimization / stress constraints / 3D modeling / FEA
© H.S. Gebremedhen et al., published by EDP Sciences, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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–15].
Most of the developments and researches on structural topology optimization are formulated and solved for minimizing compliance [15–17]. The difficulty in considering multiple loading, absence of consideration for selfweight of a structure, and stress and displacement constraints make difficult optimal layouts from this formulation to be utilized in realworld 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 twodimensional problems and those that are formulated for threedimensional problems consider a compliance minimization problem [15,20]. So far stress threedimensional 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 threedimensional optimization problem using SIMP method.
2 Problem formulation
A stressbased 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:$$\begin{array}{l}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}{\displaystyle \underset{X\mathrm{}}{\overset{}{{\displaystyle \mathrm{Min}}}}}V={\displaystyle {\displaystyle \sum}_{e=1}^{N}}{{x}_{e}}^{P}{v}_{e}\hfill \\ \text{Subjected to}:g\left({x}_{e}\right)=\frac{{\sigma}_{\mathit{\text{vm}}}}{{\sigma}_{\mathrm{\text{yield}}}}<1\hfill \\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}:KU=F\hfill \\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}0\le \rho \le 1\hfill \end{array}$$(1)
Here,
V = the volume of structure to be minimized
P_{e} = design variable
v = elemental volume
P = penalization factor
σ _{ vm } = von Mises stress induced in each element
σ _{ 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). $$\sqrt{\left[{({\sigma}_{11}{\sigma}_{22})}^{2}+{({\sigma}_{22}{\sigma}_{33})}^{2}+{({\sigma}_{11}{\sigma}_{33})}^{2}+6\left({{\sigma}_{12}}^{2}+{{\sigma}_{23}}^{2}+{{\sigma}_{13}}^{2}\right)\right]}\le {\sigma}_{\text{vm}}$$(2)
where σ _{11}, σ _{22}, and σ _{33} are the principal stresses along X, Y, and Z axes, respectively; σ _{12}, σ _{23}, and σ _{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).$$\frac{\rho {\sigma}_{ve}}{{\rho}^{p}{\sigma}_{\text{yield}}}=\xi (1\rho )$$(3)where ξ 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):$$\sigma (x)=\frac{{D}_{e}(x)\u03f5(x)}{{\rho}^{q}(x)}$$(4)where σ(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).
ϵ(x) macroscopic average strain of a material point
D_{e} (x) macroscopic elasticity tensor$${D}_{e}={\rho}^{p}{D}_{0}$$(5)
Here,
D _{0} = constitutive elasticity tensor (6):$$\left[D\right]=\frac{E}{\left(1+v\right)\left(12v\right)}\times \left[\begin{array}{c}\hfill \begin{array}{c}\hfill 1v\hfill \\ \hfill \begin{array}{c}\hfill v\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\begin{array}{c}\hfill \begin{array}{c}\hfill v\hfill \\ \hfill \begin{array}{c}\hfill 1v\hfill \\ \hfill v\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\begin{array}{c}\hfill \begin{array}{c}\hfill 0\hfill \\ \hfill \begin{array}{c}\hfill v\hfill \\ \hfill 1v\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\begin{array}{c}\hfill \begin{array}{c}\hfill 0\hfill \\ \hfill \begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill \frac{12v}{2}\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\begin{array}{c}\hfill \begin{array}{c}\hfill 0\hfill \\ \hfill \begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill \frac{12v}{2}\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\begin{array}{c}\hfill \begin{array}{c}\hfill 0\hfill \\ \hfill \begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill 0\hfill \end{array}\hfill \\ \hfill \frac{12v}{2}\hfill \end{array}\right]$$(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):$${\u03f5}_{e}={B}_{e}{u}_{e}$$(7)
B_{e} is the strain displacement matrix, which will be used to calculate stiffness matrix k _{ i }(ρ _{ e }) defined in equation (8), can be calculated as shown in equation (10) and u _{ e } is an elemental displacement vector.$${K}_{e}({\rho}_{e})={x}^{P}{k}_{e}^{o}$$(8) $${k}_{i}\left({x}_{e}\right)={\displaystyle \underset{1}{\overset{1}{{\displaystyle \int}}}}{\displaystyle \underset{1}{\overset{1}{{\displaystyle \int}}}}{\displaystyle \underset{1}{\overset{1}{{\displaystyle \int}}}}({B}^{T}DB)d\xi d\eta d\gamma $$(9)where ξ, η, and γ are the natural coordinates along X, Y, and Z axes, as shown in Figure 1, respectively.$$\left[{B}_{e}\right]=\left[\begin{array}{c}\hfill \begin{array}{ccc}\hfill \frac{\partial N1}{\partial \xi}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{\partial N1}{\partial \eta}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial N1}{\partial \gamma}\hfill \end{array}\hfill \\ \hfill \begin{array}{ccc}\hfill \frac{\partial N1}{\partial \eta}\hfill & \hfill \frac{\partial N1}{\partial \xi}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{\partial N1}{\partial \gamma}\hfill & \hfill \frac{\partial N1}{\partial \eta}\hfill \\ \hfill \frac{\partial N1}{\partial \gamma}\hfill & \hfill 0\hfill & \hfill \frac{\partial N1}{\partial \xi}\hfill \end{array}\hfill \end{array}...\begin{array}{c}\hfill \begin{array}{ccc}\hfill \frac{\partial N8}{\partial \xi}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{\partial N8}{\partial \eta}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{\partial N8}{\partial \gamma}\hfill \end{array}\hfill \\ \hfill \begin{array}{ccc}\hfill \frac{\partial N8}{\partial \eta}\hfill & \hfill \frac{\partial N8}{\partial \xi}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{\partial N8}{\partial \gamma}\hfill & \hfill \frac{\partial N8}{\partial \eta}\hfill \\ \hfill \frac{\partial N8}{\partial \gamma}\hfill & \hfill 0\hfill & \hfill \frac{\partial N8}{\partial \xi}\hfill \end{array}\hfill \end{array}\right]$$(10)
Here, N1–N8 are shape functions in natural coordinate systems defined by equation (11):$$\begin{array}{cc}\hfill {N}_{1}=\frac{\left(1\xi \right)\left(1\eta \right)\left(1\gamma \right)}{8}\hfill & \hfill {N}_{5}=\frac{\left(1\xi \right)\left(1\eta \right)\left(1+\gamma \right)}{8}\hfill \\ \hfill {N}_{2}=\frac{\left(1+\xi \right)\left(1\eta \right)\left(1\gamma \right)}{8}\hfill & \hfill {N}_{6}=\frac{\left(1+\xi \right)\left(1\eta \right)\left(1+\gamma \right)}{8}\hfill \\ \hfill {N}_{3}=\frac{\left(1+\xi \right)\left(1+\eta \right)\left(1\gamma \right)}{8}\hfill & \hfill {N}_{7}=\frac{\left(1+\xi \right)\left(1+\eta \right)\left(1+\gamma \right)}{8}\hfill \\ \hfill {N}_{4}=\frac{\left(1\xi \right)\left(1+\eta \right)\left(1\gamma \right)}{8}\hfill & \hfill {N}_{8}=\frac{\left(1\xi \right)\left(1+\eta \right)\left(1+\gamma \right)}{8}\hfill \end{array}$$(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):$$\begin{array}{l}\text{n1}=(\text{nelx+1})(\text{nely+1})(\text{elz}\text{1})\text{+}(\text{elx}\text{1})(\text{nely+1})\text{+ely}\hfill \\ \text{n2}=\text{}(\text{nelx+1})(\text{nely+1})(\text{elz}\text{1})\text{+elx}(\text{nely+1})\text{+ely}\hfill \\ \text{n3}=\text{}(\text{nelx+1})(\text{nely+1})\text{elz+elx}(\text{nely+1})\text{+ely}\hfill \\ \text{n4}=\text{}(\text{nelx+1})(\text{nely+1})\text{elz+ely}\hfill \end{array}$$(12) $$\mathrm{e}\mathrm{d}\mathrm{o}\mathrm{f}=\left[\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{ccc}\hfill 3\text{n}1+1\hfill & \hfill 3\text{n}1+2\hfill & \hfill 3\text{n}1+3\hfill \\ \hfill 3\text{n}2+1\hfill & \hfill 3\text{n}2+2\hfill & \hfill 3\text{n}2+3\hfill \\ \hfill 3\text{n}22\hfill & \hfill 3\text{n}21\hfill & \hfill 3\text{n}2\hfill \end{array}\hfill \\ \hfill \begin{array}{ccc}\hfill 3\text{n}12\hfill & \hfill 3\text{n}11\hfill & \hfill 3\text{n}1\hfill \\ \hfill 3\text{n}4+1\hfill & \hfill 3\text{n}4+2\hfill & \hfill 3\text{n}4+3\hfill \\ \hfill 3\text{n}3+1\hfill & \hfill 3\text{n}3+2\hfill & \hfill 3\text{n}3+3\hfill \end{array}\hfill \end{array}\hfill \\ \hfill \begin{array}{ccc}\hfill 3\text{n}32\hfill & \hfill 3\text{n}31\hfill & \hfill 3\text{n}3\hfill \end{array}\hfill \\ \hfill \begin{array}{ccc}\hfill 3\text{n}42\hfill & \hfill 3\text{n}41\hfill & \hfill 3\text{n}4\hfill \end{array}\hfill \end{array}\right]$$(13)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):$$KU=F$$(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):$$Ue=U\text{\hspace{0.17em}}(\mathrm{\text{edof}})$$(15)
Having the definition for strain and constitutive matrix, stress state defined in equation (4) can be expressed as equation (16):$$\sigma ({\rho}_{e})={\rho}_{e}^{Pq}{D}_{o}B{u}_{e}$$(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.$$\begin{array}{l}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\mathrm{}{\displaystyle \underset{X\mathrm{}}{\overset{}{{\displaystyle \mathrm{Min}}}}}V={\displaystyle {\displaystyle \sum}_{e=1}^{N}}\left({x}_{e}\right){v}_{e}\hfill \\ \text{Subjected to}:\text{}\frac{P{x}_{e}^{Pq}{D}_{0}B{u}_{e}}{{\sigma}_{\mathrm{y}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}}<\xi \left(1{x}_{e}\right)\hfill \\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}:KU=F\hfill \\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}0<{x}_{\text{min}}\le {x}_{e}\le 1\hfill \end{array}$$(17)
An optimality criteria method is used to solve the defined problem, by which updating scheme for design variables can be expressed as equation (18):$$\begin{array}{l}\text{if}{x}_{e}{\beta}_{e}^{\eta}\le \mathrm{max}\left({x}_{\mathrm{min}},{x}_{e}m\right)\hfill \\ {x}_{e}^{\mathrm{n}\mathrm{e}\mathrm{w}}=\mathrm{max}\left({x}_{\mathrm{min}},{x}_{e}m\right)\hfill \\ \mathrm{i}\mathrm{f}\mathrm{max}\left({x}_{\mathrm{min}},{x}_{e}m\right)<{x}_{e}{\beta}_{e}^{\eta}\le \mathrm{min}\left(1,{x}_{e}+m\right)\hfill \\ {x}_{e}^{\text{new}}={x}_{e}{\beta}_{e}^{\eta}\hfill \\ \mathrm{i}\mathrm{f}\mathrm{min}\left(1,{x}_{e}+m\right)<{x}_{e}{\beta}_{e}^{\eta}\hfill \\ {x}_{e}^{\text{new}}=\mathrm{min}\left(1,{x}_{e}+m\right)\hfill \end{array}$$(18)where m is a positive limit, which usually takes a value of 0.2,η is a numerical damping coefficient with a value of 0.5 [26,27],β is dependent on the type of problems defined as shown in equation (19):$${\beta}_{e}=\frac{\frac{\partial v}{\partial x}}{\lambda \frac{\partial g}{\partial x}}$$(19)where $\frac{\partial v}{\partial x}$ and $\frac{\partial g}{\partial x}$ are derivative of objective and constraint functions with respect to the design variable, respectively, and λ 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):$$\frac{\frac{\partial V}{\partial x}}{\frac{\partial g}{\partial x}}=\frac{{v}_{e}}{\frac{\partial \left(\frac{{\sigma}^{\mathrm{\text{vm}}}}{{\sigma}_{\mathrm{\text{yield}}}}1\right)}{\partial x}}$$(20)
The derivative of the objective V(ρ _{ e }) with respect to the design variable ρ _{ e } can be calculated as$$\frac{\partial V({\rho}_{e})}{\partial {\rho}_{e}}={\displaystyle \sum _{e\text{}=\text{}1}^{N}}\frac{\partial ({\rho}_{e}{v}_{e})}{\partial {\rho}_{e}}={\displaystyle \sum _{e\text{}=\text{}1}^{N}}{v}_{e}$$(21)
The derivative of the constraint function g(ρ _{ e }) with respect to the design variable ρ _{ e } can be calculated as$$\begin{array}{l}\frac{\partial g({\rho}_{e})}{\partial {\rho}_{e}}\Rightarrow \frac{\partial}{\partial {\rho}_{e}}(\frac{{\rho}_{e}^{Pq}{C}_{o}B{u}_{e}}{{\sigma}_{Y}}1\le {\varsigma}^{P}\varsigma )\Rightarrow \text{\hspace{0.17em}}\frac{\partial ({\rho}_{e}^{Pq})}{\partial {\rho}_{e}}\frac{{C}_{o}B{u}_{e}}{{\sigma}_{Y}}+\frac{{\rho}_{e}^{Pq}{C}_{o}B}{{\sigma}_{Y}}\frac{\partial {u}_{e}}{\partial {\rho}_{e}}\hfill \\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\hfill \end{array}$$(22)
The derivative of displacement vector can be calculated from the equilibrium equation defined in equation (14) as$$\frac{\partial K(\rho )}{\partial \rho}U(\rho )+K(\rho )\frac{\partial U(\rho )}{\partial \rho}=\frac{\partial F}{\partial \rho}=0$$(23)
which yields$$\frac{\partial U(\rho )}{\partial \rho}=K{(\rho )}^{1}\frac{\partial K(\rho )}{\partial \rho}U(\rho )$$(24)
Derivative of the stiffness matrix can be calculated from equation (8) as$$\frac{\partial {K}_{e}({\rho}_{e})}{\partial \rho}=\frac{\partial ({x}^{P}{k}_{e}^{o})}{\partial \rho}\Rightarrow p{x}^{P1}{k}_{e}^{o}$$(25)
Substituting the derivative of the stiffness matrix in equation (25), the derivative of the displacement vector becomes$$\frac{\partial U(\rho )}{\partial \rho}={({x}^{P}{k}_{e}^{o})}^{1}p{x}^{P1}{k}_{e}^{o}U(\rho )$$(26)
The derivative of the constraint function can be expressed as$$\frac{\partial g({\rho}_{e})}{\partial {\rho}_{e}}=\frac{\partial ({\rho}_{e}^{Pq})}{\partial {\rho}_{e}}\frac{{C}_{o}B{u}_{e}}{{\sigma}_{Y}}+\frac{{\rho}_{e}^{Pq}{C}_{o}B}{{\sigma}_{Y}}(p{x}^{1}U(\rho ))$$(27)
Therefore,
$${\beta}_{e}=\frac{{\displaystyle \sum _{e\text{}=\text{}1}^{N}}{v}_{e}}{\lambda \left(\frac{\partial ({\rho}_{e}^{Pq})}{\partial {\rho}_{e}}\frac{{C}_{o}B{u}_{e}}{{\sigma}_{Y}}+\frac{{\rho}_{e}^{Pq}{C}_{o}B}{{\sigma}_{Y}}(p{x}^{1}U(\rho ))\right)}$$
Fig. 1 Eightnode brick element. 
3 Numerical results
Validating the developed model benchmark problems includes simply supported beam, classical cantilever beam, and Lshape 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 threedimensional compliancebased 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 compliancebased formulation is calculated during each iteration. Compliance of a design domain for stressbased formulations is calculated during each iteration and compared with corresponding compliancebased 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 threedimensional compliancebased model developed by Liu and Tovar [15]. All benchmark problems considered in this article are simulated using MATLAB on Dell −i54 GB RAM (3.41)3.2 GHz computer.
Fig. 2 Flowchart for the developed MATLAB code to solve formulated problems. 
3.1 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 stressbased 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 STObased 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.
Fig. 3 Loading and boundary conditions for classical cantilever beam. 
Fig. 4 Optimal material distribution (a) developed model and (c) [15] and STL for optimal layout from (b) developed model and (d) [15]. 
Fig. 5 Variation of (a) maximum stress induced and (b) compliance of a cantilever beam. 
Fig. 6 Optimal material distribution (a) developed model and (b) [15] and STL for optimal layout from (c) developed model and (d) [15]. 
Fig. 7 Variation of (a) maximum stress and (b) compliance of a cantilever beam with predefined shape. 
3.2 Simply supported beam
The third case study considered is a simply supported beam under the loading and boundary conditions, as shown in 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 compliancebased topology optimization [15], respectively. As it can be seen from the figures, the optimal plots using stressbased 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 advisable to consider optimal layouts from STObased 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 threedimensional stressbased structural topology optimization.
Fig. 8 Boundary and loading conditions. 
Fig. 9 Optimal material distribution (a) developed model and (b) [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]. 
4 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 twodimensional cases.
This article presents a mathematical model for threedimensional stressbased 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.
Acknowledgments
The authors acknowledge Universiti Teknologi PETRONAS for the financial support to produce this article. The work is partially supported by Ministry of Higher Education (MOHE) Malaysia under FRGS grant FRGS Phase 1 2014.
References
 G.I.N. Rozvany, M. Zhou, T. Birker, Generalized shape optimization without homogenization. Struct. Optim. 4 , 250–252 (1992) [CrossRef] [Google Scholar]
 M.P. Bendsøe, Optimal shape design as a material distribution problem. Struct. Optim. 1 , 193–202 (1989) [CrossRef] [Google Scholar]
 K. Suzuki, N. Kikuchi, Homogenization method for shape and topology optimization. Comput. Methods Appl. Mech. Eng. 93 , 291–318 (1991) [Google Scholar]
 K. Suzuki, N. Kikuchi, A homogenization method for shape and topology optimization. Comput. Methods Appl. Mech. Eng. 93 , 291–318 (1989) [Google Scholar]
 X. Guo, et al. Stressrelated topology optimization via level set approach. Comput.Methods Appl. Mech. Eng. 200 , 3439–3452 (2011) [CrossRef] [Google Scholar]
 Q. Xia, et al. A level set solution to the stressbased structural shape and topology optimization. Comput. Struct. 90 , 55–64 (2012) [CrossRef] [Google Scholar]
 Y.M. Xie, X. Huang, Recent developments in evolutionary structural optimization (ESO) for continuum structures. IOP Conf. Ser. Mater. Sci. Eng. 10 , 012196 (2010) [CrossRef] [Google Scholar]
 Y.M. Xie, G.P. Steven, Evolutionary structural optimization (Springer, London, 1997) [Google Scholar]
 Y. Li, et al. Bidirectional evolutionary structural optimization for design of compliant mechanisms. Key Eng. Mater. 535–536 , 373–376 (2013) [CrossRef] [Google Scholar]
 O.M. Querin, G.P. Steven, Y.M. Xie, Evolutionary structural optimisation (ESO) using a bidirectional algorithm. Eng. Comput. 15 , 1031–1048 (1998) [CrossRef] [Google Scholar]
 J.D. Deaton, R.V. Grandhi, A survey of structural and multidisciplinary continuum topology optimization: post 2000. Struct. Multidiscipl. Optim. 49 , 1–38 (2014) [Google Scholar]
 G.G. Tejani, et al. Size, shape, and topology optimization of planar and space trusses using mutationbased improved metaheuristics. J. Comput. Des. Eng. 5 , 198–214 (2018) [Google Scholar]
 V.J. Savsani, G.G. Tejani, V.K. Patel, Truss topology optimization with static and dynamic constraints using modified subpopulation teachinglearningbased optimization. Eng. Optim. 48 , 1990–2006 (2016) [CrossRef] [Google Scholar]
 V.J. Savsani, et al. Modified metaheuristics using random mutation for truss topology optimization with static and dynamic constraints. J. Comput. Des. Eng. 4 , 106–130 (2017) [Google Scholar]
 K. Liu, A. Tovar, An efficient 3D topology optimization code written in MATLAB. Struct. Multidiscipl. Optim. 50 , 1175–1196 (2014) [CrossRef] [Google Scholar]
 M. Bruggi, P. Duysinx, Topology optimization for minimum weight with compliance and stress constraints. Struct. Multidiscipl. Optim. 46 , 369–384 (2012) [CrossRef] [Google Scholar]
 A. Erik, et al. Efficient topology optimization in MATLAB using 88 lines of code. Struct. Multidiscipl. Optim. 43 , 1–16 (2011) [Google Scholar]
 R. Picelli, et al. Stress minimization using the level set topology optimization, in 58th AIAA/ASCE/AHS/ASC Structures , Structural Dynamics, and Materials Conference , Grapevine, Texas, 2017 [Google Scholar]
 N. Changizi, H. Kaboodanian, M. Jalalpour, Stressbased topology optimization of frame structures under geometric uncertainty. Comput. Methods Appl. Mech. Eng. 315 , 121–140 (2017) [CrossRef] [Google Scholar]
 O. Sigmund, K. Maute, Topology optimization approaches. Struct. Multidiscipl. Optim. 48 , 1031–1055 (2013) [Google Scholar]
 H. Li, et al. A level set method for topological shape optimization of 3D structures with extrusion constraints. Comput. Methods Appl. Mech. Eng. 283 , 615–635 (2015) [CrossRef] [Google Scholar]
 P.D. Dunning, B.K. Stanford, H.A. Kim, Coupled aerostructural topology optimization using a level set method for 3D aircraft wings. Struct. Multidiscipl. Optim. 51, 1113–1132 (2015) [CrossRef] [Google Scholar]
 R. Huang, X. Huang, MATLAB implementation of 3D topology optimization using BESO, in: The 21st Australian Conference on the Mechanics of Structures and Materials. Taylor & Francis Group, London, 2011. [Google Scholar]
 P. Duysinx, O. Sigmund, New development in handling stress constraints in optimal material distribution, in: 7th AIAA /USAF/NASA/ISSMO Symposium on multidisciplinary analysis and optimization, vol. 3, St. Louis, Missouri, 1998, pp. 1501–1509. [Google Scholar]
 P. Duysinx, M.P. Bendsøe, Topology optimization of continuum structures with local stress constraints. Int. J. Numer. Methods Eng. 43 , 1453–1478 (1998) [Google Scholar]
 A. Erik, et al. Efficient topology optimization in MATLAB using 88 lines of code. Struct. Mult. Optim. 43 , 1–16 (2011) [Google Scholar]
 O. Sigmund, A 99 line topology optimization code written in MALTLAB. Struct. Multidiscipl. Optim. 21 , 120–127 (2001) [Google Scholar]
Cite this article as: Hailu Shimels Gebremedhen, Dereje Engida Woldemicahel, Fakheruldin M. Hashim, Threedimensional stressbased topology optimization using SIMP method, Int. J. Simul. Multidisci. Des. Optim. 10, A1 (2019)
All Figures
Fig. 1 Eightnode brick element. 

In the text 
Fig. 2 Flowchart for the developed MATLAB code to solve formulated problems. 

In the text 
Fig. 3 Loading and boundary conditions for classical cantilever beam. 

In the text 
Fig. 4 Optimal material distribution (a) developed model and (c) [15] and STL for optimal layout from (b) developed model and (d) [15]. 

In the text 
Fig. 5 Variation of (a) maximum stress induced and (b) compliance of a cantilever beam. 

In the text 
Fig. 6 Optimal material distribution (a) developed model and (b) [15] and STL for optimal layout from (c) developed model and (d) [15]. 

In the text 
Fig. 7 Variation of (a) maximum stress and (b) compliance of a cantilever beam with predefined shape. 

In the text 
Fig. 8 Boundary and loading conditions. 

In the text 
Fig. 9 Optimal material distribution (a) developed model and (b) [15] and STL for optimal layout from (c) developed model and (d) [15]. 

In the text 
Fig. 10 Variation of (a) maximum stress and (b) [15] and STL for optimal layout from (c) developed model and (d) [15]. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.