Basis of topological optimization

SIMP material interpolation model

Solid Isotropic Material with Penalization (SIMP) serves as the basis of topological optimization, providing an interpolation model for materials as below:
\[
E(x_i)=x_i^pE_1+(1-x_i^p)E_2
\]

\[
E(x_i=1)=E_1
\]

\[
E(x_i=x_{min})=x_{min}^pE_1+(1-x_{min}^p)E_2
\]

\(x_i^p\) is the design variable to determine how different materials are distributed or in what proportion within a single structural element. The exponent p is a penalty factor for driving the design variable approach 1 and 0. Normally p >= 3, why?

Reason 1:
Experience has proved that the punishment effect is better when p>=3.

Reason 2:
To satisfy the Hashin-Shtrikman bounds for two-phase materials. Otherwise, the equivalent stiffness generated by the SIMP interpolation model does not hold physically. (Bendsøe, M. P., 1999)

For linear elastic material, the elemental stiffness matrix can be expressed as:
\[
\mathbf{K}_i(x_i) = x_i^p\mathbf{K}_1+(1- x_i^p)\mathbf{K}_2
\]

When the Poisson ratio of two materials are equal, the elemental stiffness matrix can be written as:

Elemental sensitivity

The core of topology optimization is to determine which elements contribute more to the overall structure—retaining those with greater contributions and removing those with lesser ones. The overall structural response can be represented by compliance:

\[
C=\sum_i\mathbf{u}_i^T\mathbf{K}_i(x_i)\mathbf{u}_i
\]

To determine which element influence more on the structural compliance, the elemental sensitivity \(\alpha_i\) is constructed as below:
\[
\alpha_i =\frac{\partial C}{\partial x_i} = \mathbf{u}_i^T\frac{d\mathbf{K}_i(x_i)}{dx_i}\mathbf{u}_i = \mathbf{u}_i^T(px_i^{p-1}(\mathbf{K}_1-\mathbf{K}_2))\mathbf{u}_i=px_i^{p-1}\mathbf{u}_i^T(\mathbf{K}_1-\mathbf{K}_2)\mathbf{u}_i
\]

From another perspective, it can also be interpreted as the energy difference corresponding to the stiffness difference between the two materials.

When the Poisson ratio of two materials are equal, the elemental sensitivity can be derived as:

For solid-solid problem:

For stronger solid element:

\[
\alpha_i(x_i=1)=p(1-\frac{E_2}{E_1})\mathbf{u}_i^T\mathbf{K}_1\mathbf{u}_i
\]

For weaker solid element:

For solid-void problem:

if \(E_2\) is void, then the whole expressions degrade to the forms as:

For stronger solid element:

\[
\alpha_i(x_i=1)=p\mathbf{u}_i^T\mathbf{K}_1\mathbf{u}_i
\]

For void element:

\[
\alpha_i(x_i=x_{min})\approx 0
\]

FEA implementation

According to the general form of elemental sensitivity, \(\mathbf{u}_i^T\mathbf{K}_i(x_i)\mathbf{u}_i\) is exactly twice the elemental strain energy which can be retrieved directly from FEA software.

Since the constant coefficient does not affect the ranking of sensitivity in BESO, thus for BESO, sensitivity can be simplified to:
\[
\alpha_i=x_i^{p-1}\frac{(E_1-E_2)}{x_i^pE_1+(1-x_i^p)E_2}×strain\ energy
\]

For solid-void case, it can be further simplified to:
\[
\alpha_i=x_i^{-1}×strain\ energy
\]