# SIMP Method for Topology Optimization

Topology optimization is the most common type of structural optimization. It is used in the initial phase of the design to predict the optimal material distribution within a given initial design space of a structure and considers functional specifications and manufacturing constraints.

The most popular mathematical method for topology optimization is the Solid Isotropic Material with Penalization method (SIMP). Bendsoe and Kikuchi (1988) and Rozvany and Zhou (1992). initially proposed the SIMP method. The SIMP method predicts an optimal material distribution within a given design space, for given load cases, boundary conditions, manufacturing constraints, and performance requirements.

- ρ
_{(e)}= 1 where material is required (black) - ρ
_{(e)}= 0 where material is removed (white)

For example, the image shows an optimized material layout of a loaded
beam. The solid elements with densities ρ_{(e)} =1 are
black, whereas the void elements with ρ_{(e)} = 0 are
removed.

The introduction of a continuous relative density distribution function avoids
the binary, on-off nature of the problem. For each element, the assigned relative
density can vary between a minimum value ρ_{min} and 1,
which allows the assignment of intermediate densities for elements (characterized as
porous elements):

ρ_{min} is the minimum allowable relative density value
for empty elements that are greater than zero. This density value ensures the
numerical stability of the finite element analysis.

_{e}and the Young modulus of elasticity of the assigned isotropic material model Ε

_{0}is computed by the power law:

The penalty factor p diminishes the contribution of elements with intermediate
densities (gray elements) to the total stiffness. The penalty factor steers the
optimization solution to elements that are either solid black (ρ_{e} = 1) or void white (ρ_{e}= ρ_{min}). Numerical experiments indicate that a penalty factor
value of p = 3 is suitable.

where is the element stiffness matrix, ρ

_{min}is the minimum relative density, ρ

_{e}is the element relative density, p is the penalty factor, and N is the number of elements in the design domain.

For example, for an element with an assigned relative density ρ_{e} = 0.5, penalty factor = 3, and ρ_{min} = 0.001, the global stiffness matrix is scaled by a factor of (0.001
+ (1 -0.001)* 0.5 ^3) = 0.12587.

## Objective Function: Maximizing Stiffness

A popular optimization objective is to maximize the overall stiffness of a structure, or minimize its compliance under a given amount of mass removal.

Compliance is a measure of the overall flexibility or softness of a structure, and is the reciprocal of stiffness. The global compliance is equal to the sum of the element elastic or strain energies. Minimizing the global compliance, C, is equivalent to maximizing the global stiffness. The optimization algorithm, through an iterative process, seeks to resolve the element densities (which are the optimization design variables) that minimize the global compliance of the structure.

[u_{e}] is the nodal displacement vector of element e, [K_{e}] is the stiffness
of element e, and vector {ρ} contains the elements' relative densities
**ρ**_{e}.

v

_{e}is the element volume, and M

_{target}is the target mass of the optimization.

[K{

**ρ**}] is the global stiffness matrix modulated by the vector of relative densities, {u} is the displacement vector, and {F} is the external force vector.

The formula above contains design response constraints such as limits on stresses, displacements, eigenfrequencies, etc.

## Sensitivity Analysis

During each iteration, the optimization algorithm performs a sensitivity analysis to evaluate the impact the variation of material densities has on the objective function to maximize stiffness.

Mathematically, the sensitivity analysis is expressed as the derivative of the objective function with respect to the material densities:

During a sensitivity analysis, elements weighted with low material density factors eventually lose their structural importance and are eliminated during further iterations.

If you calculate the sensitivity for each element independently and don't consider the connectivity between elements, this can lead to material discontinuity and to volumes becoming unconnected to the main geometry. This is known as checkerboard effect. To reduce the checkerboard effect, a filtering scheme applies an element influence radius and averages the sensitivities of each element inside its influence region.

The optimization iterations continue until the variations of the objective function converge and iterations reach their converge criteria.