Stress-Constrained Topology Optimization for Commercial Software: A Python Implementation for ABAQUS^®

Abstract

Topology optimization has evidenced its capacity to provide new optimal designs in many different disciplines. However, most novel methods are difficult to apply in commercial software, limiting their use in the academic field and hindering their application in the industry. This article presents a new open methodology for solving geometrically complex non-self-adjoint topology optimization problems, including stress-constrained and stress minimization formulations, using validated FEM commercial software. The methodology was validated by comparing the sensitivity analysis with the results obtained through finite differences and solving two benchmark problems with the following optimizers: Optimality Criteria, Method of Moving Asymptotes, Sequential Least-Squares Quadratic Programming (SLSQP), and Trust-constr optimization algorithms. The SLSQP and Trust-constr optimization algorithms obtained better results in stress-minimization problem statements than the methodology available in ABAQUS^®. A Python implementation of this methodology is proposed, working in conjunction with the commercial software ABAQUS^® 2023 to allow a straightforward application to new problems while benefiting from a graphic user interface and validated finite element solver.

1. Introduction

Originally described by Bendsøe and Kikuchi in 1988 [1], topology optimization is one of the three sub-fields of structural optimization, amongst size and shape optimization. Usually applied in the early stage of structural design, its purpose is to find the optimal distribution of material when certain load conditions are applied [2,3]. The design region is usually represented with density based or level-set functions, which are optimized according to a certain objective function and constraints, defining which regions should have material and which should not [4].

The popularity of topology optimization has led to the publishing of several educational articles [5]. These include the MATLAB codes written by Ole Sigmund [6], later revisited by Andreassen et al. [7], and an equivalent implementation for 3D problems by Liu [8], all of them using the Solid Isotropic Material with Penalization (SIMP) method [9,10] to define the material properties and the Optimality Criteria (OCs) [11] to determine the optimal design variables. Other publications include the MATLAB implementation of the level-set method [12,13] by Challis [14] and the implementations of the Bi-directional Evolutionary Structural Optimization (BESO) method [15] in MATLAB and Python by Huang and Xie [16] and Zuo and Xie [17], respectively.

Table 1. Literature review summary table, including scientific articles referred to in the ABAQUS^® and ANSYS^® documentation, sorted by scope. The black dots indicate an affirmative answer to the questions posed in each column.

Reference Source	Scope of the Reference	Author	Algorithm	Provides Code?	Consider Non-Material Constraint?	Suitable for Complex Geometries?	Based on Commercial Solver?
Literature	Review article or complementary method	Sigmund [18]	-
		Zhu et al. [19]	-
		Sigmund and Maute [20]	OC	•
		Wang et al. [5]	-
		Harzheim and Graf [21]	-
	Compliance minimization	Sigmund [6]	OC	•
		Andreassen et al. [7]	OC	•
		Challis [14]	Level-set method	•
		Zuo and Xie [17]	OC	•			•
		Suresh [22]	Pareto-optimal tracing	•
		Talischi et al. [23]	OC	•
		Smith and Norato [24]	MMA and fmincon	•
		Sanders et al. [25]	Modified OC	•
		Otomori et al. [26]	Level-set method	•
		Gao et al. [27]	OC	•
		Xia and Breitkopf [28]	OC	•
		Liu and Tovar [8]	OC, SLSQP, and MMA	•
		Aage et al. [29]	MMA	•		•
		Liang and Cheng [30]	OC	•
		Picelli et al. [31]	OC	•
		Ferrari and Sigmund [32]	OC	•
		Du et al. [33]	Level-set method	•
		Chen et al. [34]	MMA	•			•
		Sotiropoulos et al. [35]	OC, MMA
	Stress-constrained compliance minimization	Londoño and Paulino [36]	MMA	•	•
		Yang et al. [37]	MMA		•
		Bruggi and Venini [38]	MMA		•
		Lee et al. [39]	SQP		•
		Cai and Zhang [40]	Level-set method		•
		Troya and Tortorelli [41]	MMA		•
		Suresh and Takalloozadeh [42]	Level-set method		•	•
		Chu et al. [43]	OC		•	•
		Oh et al. [44]	Gradient-based phase field		•
		Luo et al. [45]	MMA		•
		Biyikli and To [46]	PTO	•	•
		Amir [47]	MMA	•	•
		París et al. [48]	SLP and steepest descent		•
		Holmberg et al. [49]	MMA		•
		Norato et al. [50]	MMA		•
		Burger and Stainko [51]	Interior point method		•
		De Leon et al. [52]	MMA		•
		Holmberg et al. [53]	MMA		•
		Pastore et al. [54]	SLP, MMA, and OC		•
		París et al. [55]	SLP and steepest descent		•
		Deng and Suresh [56]	Level-set method		•
		Senhora et al. [57]	MMA		•
		Saadlaoui et al. [58]	Commercial software (black box)		•	•	•
		Holmberg et al. [59]	MMA		•
		Amir and Lazarov [60]	MMA		•
		Granlund et al. [61]	MMA		•
		Amstutz et al. [62]	Level-set method		•
		Cheng et al. [63]	MMA		•
		Deng et al. [64]	MMA	•	•
		Mirzendehdel and Suresh [65]	Level-set method		•	•
	Stress- and compliance-constrained mass minimization	Collet et al. [66]	MMA		•
		Bruggi and Duysinx [67]	MMA		•
	Stress-constrained mass minimization	Londoño et al. [68]	MMA		•
		París et al. [69]	SLP		•
		Long et al. [70]	SQP and MMA		•
		Duysinx [71]	MMA		•
		Silva et al. [72]	Level-set method		•	•	•
	Stress-constrained compliance maximization	Miyajima et al. [73]	Level-set method		•
		Emmendoerfer et al. [74]	Level-set method		•
ABAQUS® documentation	Review article or complementary method	BendsØe and Sigmund [75]	OC, MMA, SLP, and level-set method	•	•
	Compliance minimization	Bakhtiary et al. [76]	OC			•	•
ANSYS® documentation	Review article or complementary method	Mlejnek [77]	-
		Stolpe and Svanberg [78]	-
		Pedersen and Allinger [79]	-
		BendsØe and Sigmund [80]	-
		Clausen and Pedersen [81]	-
		Svanberg [82]	-
	Compliance minimization	Kabus and Pedersen [83]	Commercial software (black box)		•	•	•
	Minimization of maximum pressure	Søndergaard and Pedersen [84]	Commercial software (black box)		•	•	•
	Eigenvector minimization	Hansen [85]	MMA
		Olhoff and Du [86]	MMA			•	•
	Compliance and eigenfrequency minimization	Jog [87]	OC
		Sigmund and Jensen [88]	MMA		•
	Maximization of the magnitude of steady-state vibrations	Tcherniak [89]	MMA		•

compiles the publications on the development and implementation of topology optimization methods, sorted by the problem statement addressed, highlighting:

• The optimization algorithm used;
• If the publication provides the code implementation;
• If the method considers more constraints than just a mass or volume limit;
• If the method is suitable for complex geometries, allowing an easy application to case studies with intricate design spaces or the interaction with CAD models;
• If the method is based on commercial solvers.

Table 1 includes several publications with open access implementations of topology optimization methods suitable for compliance-minimization problems [6,7,14,17,22,23,24,25,26,27,28,29,30,31,32,34]. However, only four publications provide implementations of topology optimization methods compatible with stress-constrained compliance-minimization problems [36,46,47,64]. Regarding these three publications, Biyikli et al. [46] proposed a non-sensitivity method and Amir [47] proposed a method based on inexact sensitivities, leaving the MATLAB implementations proposed by Londoño and Paulino [36] and by Deng et al. [64] as the only publications providing a non-approximated solution to the stress-constrained compliance-minimization problem. On top of the reduced number of publications providing open access solutions with non-approximated methods to this problem statement, there are additional points of improvement that could be addressed to promote the research and application of topology optimization methods in the industry: including the compatibility with complex design spaces and geometries, as well as the interaction with commercially available and validated solvers. Finally, it is also relevant to consider that the topology optimization methods available in commercial software, such as ABAQUS^® 2023 and ANSYS^® 2023, operate on an opaque black-box system, limiting the access to relevant information and hindering their extension to a new feature of interest, such as reliability-based topology optimization [90,91,92] and robust topology optimization considering uncertainties [93].

To address this issue, the present article contributes with a methodology capable of solving non-self-adjoint problems in commercial software, including stress-constrained compliance minimization and stress-minimization problem statements, defined in a complex geometrical space. This methodology has been validated, evaluating the sensitivity analysis with finite differences and solving two benchmark problems with the following optimizers: Optimality Criteria, Method of Moving Asymptotes, Sequential Least-Squares Quadratic Programming, and Trust-constr optimization algorithms. In compliance-minimization problems and stress-constrained compliance-minimization problems, the methodology herein implemented achieved a performance equal to ABAQUS^® in terms of the objective function. In stress-minimization problems, the methodology herein implemented allowed a stress reduction 35% superior to the one obtained with ABAQUS^®.

To allow straightforward use even in an industrial field, a Python implementation is provided as Supplementary Material. This implementation works in conjunction with ABAQUS^®, allowing a direct application to new problems, providing a graphic user interface and access to a commercially validated finite-element solver. The implementation can also be downloaded in digital format from the following Dataset [94] or repository: https://github.com/pnfernandes/Python-Code-for-Stress-Constrained-Topology-Optimization-in-ABAQUS.

The following sections of the article are organized as follows. Section 2 and Section 3 formulate the topology optimization problem statements and deduce the necessary function derivatives in both continuous and discrete formulations, respectively. Section 4 describes the filter technique to void numerical instabilities. Section 5 overviews the optimization algorithms available in the provided code. Section 6 summarizes the most-relevant aspects of the Python code implementation, while in Section 7, the benchmark problems considered are presented. Section 8 demonstrates the suitability of the code implementation to solve the problems proposed. The results of the topology optimization processes are presented in Section 9. Finally, Section 10 summarizes the main conclusions of the present work.

2. Continuous Formulation of Topology Optimization Problem Statements and Sensitivities

In this section, the theoretical part of the topology optimization problem statement and the sensitivity analysis are defined in their continuous versions. Its practical counterpart, the discrete formulation, is presented in Section 3, which may allow an easier understanding of both the theoretical and implementation fundamentals of the topology optimization method.

2.1. Topology Optimization Problem

Consider the formulation of a generic topology optimization problem, defined as finding $\rho \in L^{\infty }\left(\mathsf{\Omega }\right)$ such that:

$$\underset{\rho }{min}:J(\rho ,u\left(\rho \right))$$

(1)

subject to:

$$\begin{array}{c}\hfill {\rho }_{min}\le \rho \le 1\end{array}$$

(2)

$$\begin{array}{c}\hfill {\int }_{\mathsf{\Omega }}\rho dx\le V^{*}\end{array}$$

(3)

where $\mathsf{\Omega }$ is the domain of the topology optimization problem, $L^{\infty }\left(\mathsf{\Omega }\right)$ is the space of bounded functions, and the density $\rho$ represents the design variables, which can vary within the interval $[{\rho }_{min},1]$. J represents a generic objective or cost function; $\int \rho dx$ represents a generic volume constraint function, where $V^{*}$ is its maximum allowable value; $u\left(\rho \right)$ are the displacement solutions of the standard equilibrium equation presented as finding $u\in H_0^1\left(\mathsf{\Omega }\right)$ such that:

$$\begin{array}{cc}\hfill a(\rho ,u,v)& =l(\rho ,v),\forall v\in H_0^1\left(\mathsf{\Omega }\right)\hfill \end{array}$$

(4)

where $a(\rho ,u,v)$ is the constitutive function of the equilibrium equation, $l(\rho ,v)$ are the external forces, and $H_0^1\left(\mathsf{\Omega }\right)$ is the space of functions with square integrable derivatives and homogeneous values on the boundary and domain $\mathsf{\Omega }\in R^D$, with $D\in [2,3]$. The bilinear form is:

$$\begin{array}{cc}\hfill a(\rho ,u,v)& ={\int }_{\mathsf{\Omega }}{▽}^{s}u\mathbb{C}\left(\rho \right){▽}^{s}vdx\hfill \end{array}$$

(5)

with the external forces $l(\rho ,v)$ written as:

$$\begin{array}{cc}\hfill l(\rho ,v)& ={\int }_{\mathsf{\Omega }}f\left(\rho \right)vdx+{\int }_{\partial \mathsf{\Omega }}t·nvds.\hfill \end{array}$$

(6)

Here, ${▽}^s$ represents a symmetric gradient and $\mathbb{C}\left(\rho \right)$ is the constitutive tensor, and thus, the strains $\xi \left(u\right)={▽}^{s}u$ are in accordance with linear elasticity. v is the corresponding test function; $f\left(\rho \right)$ is a generic volumetric load function; t indicates the applied boundary forces, if they exist; n is the normal direction pointing outwards of the boundary $\partial \mathsf{\Omega }$ of v. It is important to note that this material constraint can be applied to either the volume or mass of the structure, following the same expressions. The volume constraint is preferred and used in this article to conform with the most-used term in the literature.

When considering a maximum stress constraint, the problem statement is also subject to:

$$\begin{array}{c}\hfill {\int }_{\mathsf{\Omega }}{\sigma }_a^{VM}(\rho ,u\left(\rho \right))dx\le {\sigma }^{*}\end{array}$$

(7)

where ${\int }_{\mathsf{\Omega }}{\sigma }_a^{VM}(\rho ,u\left(\rho \right))dx$ represents the stress constraint with a maximum allowable value ${\sigma }^{*}$. In particular, ${\sigma }_a^{VM}$ is the penalized von Mises stress, as described in the following subsection.

2.2. Regularization and Penalization

Regarding this particular work and implementation, it is important to note two adopted considerations. The first one is using two penalization factors applied to the material stiffness and stress. These factors aim at making intermediate design solutions uncompetitive and, in turn, promoting black-and-white solutions.

The material stiffness is penalized by the SIMP penalization parameter $P=3.0$, in accordance with [80], when using continuous design variables (Equation (8)). In this implementation, this factor is also applied to the other material properties, as:

$$\mathbb{C}\left(\rho \right)={\rho }^P{\mathbb{C}}_0$$

(8)

where ${\mathbb{C}}_0$ is the material stiffness of a fully solid element. Note that ${\mathbb{C}}^{\prime }\left(\rho \right)=P{\rho }^{P-1}{\mathbb{C}}_0$ for this type of stiffness penalization.

The stress is also penalized according to a factor equal to ${\rho }^{\beta }$, with $\beta =\frac{1}{2}$, as adopted in [53], following the initial proposal by Bruggi et al. [95] with the exponent suggested in [96]. This penalization leads to a non-physical stress for intermediate design densities, but not for black-and-white solutions, and tends towards 0 when the design density decreases, avoiding singularity problems [53]. Thus, ${\sigma }_a$ is defined as the amplified stress, described by the following expression:

$${\sigma }_a\left(\rho \right)={\rho }^{\beta }\widehat{\sigma }\left(\rho \right)$$

(9)

with:

$$\widehat{\sigma }\left(\rho \right)={\mathbb{C}}_0{▽}^{s}u\left(\rho \right)$$

(10)

where $\widehat{\sigma }\left(\rho \right)$ is the stress vector, written in Voigt notation. Therefore, the von Mises amplified stress norm is ${\sigma }_a^{VM}={\left({\sigma }_{a}M{\sigma }_a\right)}^{\frac{1}{2}}$, with M being the von Mises matrix operator.

The second consideration is the use of a regularization approach, where removed or void elements will maintain a minimum design density, ${\rho }_{min}=0.01$ by default, which contrasts with the complete element removal. In the literature, these regularization techniques may be referred to as “soft-kill” and “hard-kill” approaches [17]. The regularization approach was introduced in the code implementation to avoid the stiffness matrix becoming singular, as well as allowing a constant mesh in the finite-element model during the topology optimization loop [17,53].

2.3. Sensitivity Analysis

The derivative of the cost function in the direction $\tilde{\rho }\in L^{\infty }\left(\mathsf{\Omega }\right)$ is then defined as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill DJ(\rho & ,u\left(\rho \right))\tilde{\rho }=D_{\rho }J(\rho ,u\left(\rho \right))\tilde{\rho }+D_{u}J(\rho ,u\left(\rho \right))D_{\rho }u\left(\rho \right)\tilde{\rho },\forall \tilde{\rho }\in L^{\infty }\left(\mathsf{\Omega }\right)\hfill \end{array}\end{array}$$

(11)

Taking the derivative in the equilibrium equation, in order to find $D_{u}J(\rho ,u\left(\rho \right))D_{\rho }u\left(\rho \right)\tilde{\rho }$, leads to:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill [D_{\rho }a& (\rho ,u\left(\rho \right),v)-D_{\rho }l(\rho ,v)]\tilde{\rho }+D_{u}a(\rho ,u\left(\rho \right),v)D_{\rho }u\left(\rho \right)\tilde{\rho }=0\hfill \\ & \forall v\in H_0^1\left(\mathsf{\Omega }\right),\forall \tilde{\rho }\in L^{\infty }\left(\mathsf{\Omega }\right).\hfill \end{array}\end{array}$$

(12)

Since $a(\rho ,u(\rho ),v)$ is linear, it can be rewritten as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill D_{u}a(\rho ,u\left(\rho \right),v)D_{\rho }u\left(\rho \right)\tilde{\rho }=a(\rho ,D_{\rho }u\left(\rho \right)\tilde{\rho },v);\end{array}\end{array}$$

(13)

therefore, Equation (12) can be rewritten as follows:

$$\begin{array}{c}\begin{array}{c}-a(\rho ,D_{\rho }u\left(\rho \right)\tilde{\rho },v)=[D_{\rho }a(\rho ,u\left(\rho \right),v)-D_{\rho }l(\rho ,v)]\tilde{\rho },\forall v,\tilde{\rho }.\hfill \end{array}\end{array}$$

(14)

Solving Equation (14) for all values of $\tilde{\rho }$ would be too expensive. For this reason, the use of the adjoint method is preferred. To do so, the adjoint variable $\lambda$ is defined, the solution of:

$$\begin{array}{c}\hfill \begin{array}{cc}& a(\rho ,\lambda ,w)=-D_{u}J(\rho ,u\left(\rho \right))w\forall w\in H_0^1\left(\mathsf{\Omega }\right)\iff \hfill \\ & \iff -a(\rho ,D_{\rho }u\left(\rho \right)\tilde{\rho },\lambda )=D_{u}J(\rho ,u\left(\rho \right))D_{\rho }u\left(\rho \right)\tilde{\rho },\forall \tilde{\rho }.\hfill \end{array}\end{array}$$

(15)

Then, taking $w=D_{\rho }u\left(\rho \right)\tilde{\rho }$ in Equation (15):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_{u}J& (\rho ,u\left(\rho \right))D_{\rho }u\left(\rho \right)\tilde{\rho }=-a(\rho ,\lambda ,D_{\rho }u\left(\rho \right)\tilde{\rho })=\hfill \\ & =-a(\rho ,D_{\rho }u\left(\rho \right)\tilde{\rho },\lambda )=\hfill \\ & =[D_{\rho }a(\rho ,u\left(\rho \right),\lambda )-D_{\rho }l(\rho ,v)]\tilde{\rho }\hfill \end{array}\end{array}$$

(16)

where the self-adjoint property of $a(\rho ,\lambda ,w)$ and Equation (12) are applied.

The generic optimization process, which can be defined in four main steps, is proposed as follows:

• Find $u\left(\rho \right)\in H_0^1\left(\mathsf{\Omega }\right)$, the solution of: $a(\rho ,u,v)=l(\rho ,v),\forall v\in H_0^1\left(\mathsf{\Omega }\right)$;
• Find $\lambda$, the solution of: $a(\rho ,\lambda ,w)=-D_{u}J(\rho ,u\left(\rho \right))w,\forall w$;
• Compute the derivative as $DJ(\rho ,u\left(\rho \right))\tilde{\rho }=[D_{\rho }J(\rho ,u\left(\rho \right))+a_{\rho }(\rho ,u\left(\rho \right),\lambda )-l_{\rho }(\rho ,\lambda )]\tilde{\rho }$;
• Update the design variables using the gradient information.

The first two steps correspond to solving the state and adjoint problems, respectively. The third step is the gradient calculation, which is defined in terms of the state and adjoint variables. The fourth and final step consists of using the gradient to determine the next value of the design variables, which can be performed by any suitable optimization algorithm (such as the algorithms described in Section 5).

2.4. Compliance Functional

In the particular problem statement of compliance minimization, the objective function can be rewritten as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill J(\rho ,u\left(\rho \right))={\int }_{\mathsf{\Omega }}fudx\end{array}\end{array}$$

(17)

and its derivative in the direction $w\in H_0^1\left(\mathsf{\Omega }\right)$ is then:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill D_{u}J(\rho ,u\left(\rho \right))w={\int }_{\mathsf{\Omega }}fwdx=l\left(w\right).\end{array}\end{array}$$

(18)

Notice that, in the particular case of compliance-minimization problem statements, it is commonly assumed in the literature that the external loads do not depend on $\rho$. Therefore:

$$\begin{array}{cc}\hfill a(\rho ,u,v)& =l\left(v\right).\hfill \end{array}$$

(19)

In this particular case, the adjoint variable leads to the same expression shown in Equation (15), and introducing Equation (4), it is implied that:

$$\begin{array}{c}\hfill \begin{array}{cc}& a(\rho ,\lambda ,w)=-D_{u}J(\rho ,u\left(\rho \right))w,\forall w\hfill \\ & \iff -a(\rho ,u,w)=-D_{u}J(\rho ,u\left(\rho \right))=-l\left(w\right).\hfill \end{array}\end{array}$$

(20)

Note that Equation (20), which describes the adjoint problem of the compliance-minimization problem statement, returns the definition of the state problem, setting them equal to each other. Functionals that are “self-adjoint” lead to clear computational benefits since Equations (4) and (15) can be solved with the same computation, a single FEA in the case of the code provided.

Since $l\left(v\right)$ and $l\left(w\right)$ do not depend on $\rho$, the terms $D_{\rho }J(\rho ,u\left(\rho \right))\tilde{\rho }=D_{\rho }l\left(w\right)\tilde{\rho }=D_{\rho }l\left(v\right)\tilde{\rho }=0$, and thus, the derivative is just:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill DJ(\rho ,u\left(\rho \right))\tilde{\rho }=-D_{\rho }a(\rho ,u\left(\rho \right),u\left(\rho \right))\tilde{\rho }.\end{array}\end{array}$$

(21)

Following the definition of the compliance, the derivative of the bilinear form is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill D_{\rho }a(\rho ,u\left(\rho \right),v)={\int }_{\mathsf{\Omega }}{▽}^{s}v{\mathbb{C}}^{\prime }\left(\rho \right){▽}^{s}u\tilde{\rho }dx\end{array}\end{array}$$

(22)

Finally, considering the stiffness penalization used, Equation (22) can be introduced into (21), obtaining:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill DJ(\rho ,u\left(\rho \right))\tilde{\rho }& ={\int }_{\mathsf{\Omega }}{▽}^{s}u\left(P{\rho }^{P-1}{\mathbb{C}}_0\right){▽}^{s}u\tilde{\rho }dx={\int }_{\mathsf{\Omega }}g\tilde{\rho }dx\hfill \end{array}\end{array}$$

(23)

where g is the gradient of the compliance objective function, thus obtained as $g={▽}^{s}u$$\left(P{\rho }^{P-1}{\mathbb{C}}_0\right){▽}^{s}u$. For this particular case, the first three optimization steps can be rewritten as:

• Find $u\left(\rho \right)\in H_0^1\left(\mathsf{\Omega }\right)$, the solution of: $a(\rho ,u\left(\rho \right),v)=l\left(v\right),\forall v\in H_0^1\left(\mathsf{\Omega }\right)$;
• Take $\lambda =-u$, since $a(\rho ,\lambda ,w)=-D_{u}J(\rho ,u\left(\rho \right))w=-l\left(v\right),\forall w$ is exactly the same problem as the first step;
• Compute: $DJ(\rho ,u\left(\rho \right))\tilde{\rho }=-a_{\rho }(\rho ,u\left(\rho \right),u\left(\rho \right))\tilde{\rho }={\int }_{\mathsf{\Omega }}{▽}^{s}u\left(P{\rho }^{P-1}{\mathbb{C}}_0\right){▽}^{s}u\tilde{\rho }dx$.

2.5. Stress Functional

In the stress-minimization or stress-constrained compliance-minimization problem statements, the derivative of the stress norm functional should also be taken. Here, the maximum function is approximated by means of a modified p-norm function. This approximation is necessary to provide a derivable function that approximates the maximum function, which is non-differentiable. Adopting the modified p-norm approximation proposed in [53], the maximum function represented in Equation (7) can be redefined as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill J(\rho ,u\left(\rho \right))={\left({\int }_{\mathsf{\Omega }}{\left({\sigma }_a^{VM}(\rho ,u\left(\rho \right))\right)}^{q}dx\right)}^{\frac{1}{q}}\end{array}\end{array}$$

(24)

where q is the exponential factor, and the von Mises amplified stress norm is ${\sigma }_a^{VM}={\left({\sigma }_{a}M{\sigma }_a\right)}^{\frac{1}{2}}$, with M being the von Mises matrix operator. Note that, although unconventional, this manuscript adopted q as the exponential factor of the p-norm approximation to improve the readability of the equation and to avoid confusion with the design density symbol $\rho$ and the SIMP exponential factor P. Furthermore, it is relevant to highlight that, by adopting the use of this approximation, a global approach is selected, instead of local enforcement at each integration point. This option was preferred for providing better scaling with the mesh refinement, reducing the computational cost, and for its simplicity. However, the use of augmented Lagrangian methods could be seen as an alternative [57].

The term $D_{\rho }J(\rho ,u\left(\rho \right))\left(\tilde{\rho }\right)$ is defined as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_{\rho }& J(\rho ,u\left(\rho \right))\left(\tilde{\rho }\right)=\hfill \\ & ={\left({\int }_{\mathsf{\Omega }}{\left({\sigma }_a^{VM}(\rho ,u\left(\rho \right))\right)}^{q}dx\right)}^{\frac{1}{q}-1}\times q{\int }_{\mathsf{\Omega }}{\sigma }_a^{VM}{(\rho ,u\left(\rho \right))}^{q-1}{D}_{{\sigma }_a}{\sigma }_a^{VM}\left(D_{\rho }{\sigma }_a\left(\tilde{\rho }\right)\right)dx\hfill \end{array}\end{array}$$

(25)

where:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill D_{\rho }{\sigma }_a\left(\tilde{\rho }\right)=\beta {\rho }^{\beta -1}{\mathbb{C}}_0{▽}^{s}u\tilde{\rho }.\end{array}\end{array}$$

(26)

The derivative in the direction w, which allows us to define the adjoint problem, is:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_{u}J& (\rho ,u(\rho \left)\right)\left(w\right)=\hfill \\ & ={\left({\int }_{\mathsf{\Omega }}{\left({\sigma }_a^{VM}(\rho ,u\left(\rho \right))\right)}^{q}dx\right)}^{\frac{1}{q}-1}\times {\int }_{\mathsf{\Omega }}{\sigma }_a^{VM}{\left(\rho ,u\left(\rho \right)\right)}^{q-1}{D}_{{\sigma }_a}{\sigma }_a^{VM}\left(D_u{\sigma }_a\left(w\right)\right)dx\hfill \end{array}\end{array}$$

(27)

where:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_{{\sigma }_a}{\sigma }_a^{VM}\left(D_u{\sigma }_a\left(w\right)\right)& ={\left({\sigma }_{a}M{\sigma }_a\right)}^{-\frac{1}{2}}{\sigma }_{a}M{D}_u{\sigma }_a\left(w\right)\hfill \end{array}\end{array}$$

(28)

with:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill D_u{\sigma }_a\left(w\right)& ={\rho }^{\beta }{\mathbb{C}}_0{▽}^{s}w.\hfill \end{array}\end{array}$$

(29)

Finally, before determining the value of $DJ(\rho ,u\left(\rho \right))\tilde{\rho }$, it is required to find $D_{\rho }a(\rho ,u\left(\rho \right),\lambda )$ as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill D_{\rho }a(\rho ,u\left(\rho \right),\lambda )\tilde{\rho }={\int }_{\mathsf{\Omega }}{▽}^s\lambda {\mathbb{C}}^{\prime }\left(\rho \right){▽}^{s}u\tilde{\rho }.\end{array}\end{array}$$

(30)

As stated previously, note that ${▽}^s\lambda$ and ${▽}^{s}u$ represent the strains of the adjoint and state problems, respectively. Therefore, in the provided code implementation, these terms are obtained directly from the ABAQUS^® FEA. Also, note that, for the stress functional, $D_{\rho }l\left(\lambda \right)=0$. Finally, the first three optimization steps can then be rewritten as:

• Find u, the solution of $a(\rho ,u(\rho \left)v\right)=l(\rho ,v)$;
• Using Equation (27), find $\lambda$ such that $a(\rho ,\lambda ,w)=-D_{u}J(\rho ,u\left(\rho \right))w,\forall w$. Note that this functional is not self-adjoint, leading to a different procedure than the one shown in the previous section;
• Compute: $DJ(\rho ,u\left(\rho \right))\tilde{\rho }=[D_{\rho }J(\rho ,u\left(\rho \right))+D_{\rho }a(\rho ,u\left(\rho \right),\lambda )]\tilde{\rho }$.

3. Discrete Formulation of Topology Optimization Problem Statements and Sensitivities

In this section, the topology optimization problem statement and sensitivity analysis defined in Section 2 are now presented in their discretized versions. This information is included here to bridge the gap between the continuous framework shown in Section 2 and the implementation used in the code provided with this work, allowing an easier understanding of both the theoretical and implementation fundamentals of a topology optimization method.

3.1. Topology Optimization Problem

The discrete version of Problem (1) results in:

$$\underset{\rho }{min}:J(\rho ,u\left(\rho \right))$$

(31)

subject to:

$${\rho }_{min}\le {\rho }_e\le 1\forall e=1,\dots ,N$$

(32)

$$V\left(\rho \right)=\sum _{\rho }{\rho }_e{v}_e\le V^{*}$$

(33)

where, as previously defined, $J\left(\rho \right)$ is the cost function, the density $\rho$ represents the design variables, which can vary within the interval $[{\rho }_{min},1]$, $V\left(\rho \right)$ is the total volume, with $v_e$ representing the volume of element e when the design density $\rho =1$, $V^{*}$ is the maximum value of the volume constraint, and $u\left(\rho \right)$ is the solution displacement vector of the state problem:

$$F=K\left(\rho \right)u.$$

(34)

The stiffness matrix $K\left(\rho \right)$ and force vector F are defined as, respectively:

$$K\left(\rho \right)={\int }_{\mathsf{\Omega }}{B}_a^T\mathbb{C}{B}_{a}dx$$

(35)

$$F={\int }_{\mathsf{\Omega }}Nfdx+{\int }_{\partial \mathsf{\Omega }}Nt·nds$$

(36)

where $B_a$ is the strain–displacement matrix in the evaluation point a, N a linear shape function, and f a generic load function, while t indicates the applied boundary forces, if they exist, and n is the normal direction pointing outwards of the boundary $\partial \mathsf{\Omega }$.

When considering a maximum stress constraint, the problem statement is also subject to:

$${\sigma }^{PN}\left(\rho \right)\le {\sigma }^{*}$$

(37)

where ${\sigma }^{PN}\left(\rho \right)$ represents the maximum stress and ${\sigma }^{*}$ its the maximum allowable value.

Note that the penalization factors and regularization process described in Section 2.2 are not changed.

3.2. Sensitivity Analysis

The gradient of the cost function can be defined as:

$${▽}_{\rho }J=\frac{\partial J}{\partial \rho }+\frac{\partial J}{\partial u}\frac{\partial u}{\partial \rho }.$$

(38)

In order to find $\frac{\partial J}{\partial u}$, it is possible to reorganize the state equation, multiply it for a generic vector v, and derive the expression, leading to:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill v^T\left[\frac{\partial K\left(\rho \right)}{\partial \rho }u\left(\rho \right)-\frac{\partial F\left(\rho \right)}{\partial \rho }+K\left(\rho \right)\frac{\partial u\left(\rho \right)}{\partial \rho }\right]=0;\end{array}\end{array}$$

(39)

therefore:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill -{\left(\frac{\partial u\left(\rho \right)}{\partial \rho }\right)}^T{K}^T\left(\rho \right)v=v^t\left[\frac{\partial K\left(\rho \right)}{\partial \rho }u-\frac{\partial F\left(\rho \right)}{\partial \rho }\right].\end{array}\end{array}$$

(40)

As stated in Section 2, solving Equation (40) for all values of $\rho$ would be too expensive, motivating the use of the adjoint method. Defining the adjoint variable as $\lambda =v$ such that $K^T\left(\rho \right)\lambda =-\frac{\partial J}{\partial u}$ leads to:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill -{\left(\frac{\partial u\left(\rho \right)}{\partial \rho }\right)}^T& K^T\left(\rho \right)\lambda ={\left(\frac{\partial u}{\partial \rho }\right)}^T\frac{\partial J}{\partial u}={\lambda }^T\left[\frac{\partial K\left(\rho \right)}{\partial \rho }u-\frac{\partial F\left(\rho \right)}{\partial \rho }\right].\hfill \end{array}\end{array}$$

(41)

Thus, the gradient can finally be rewritten as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {▽}_{\rho }J& =\frac{\partial J}{\partial \rho }+\frac{\partial J}{\partial u}\frac{\partial u}{\partial \rho }=\frac{\partial J}{\partial \rho }+{\lambda }^T\left[\frac{\partial F\left(\rho \right)}{\partial \rho }-\frac{\partial K\left(\rho \right)}{\partial \rho }u\right].\hfill \end{array}\end{array}$$

(42)

With this information, the equivalent three-step process described in Section 2.3 can be rewritten as:

• Find u such that: $K\left(\rho \right)u=F\left(\rho \right)$;
• Find $\lambda$ such that: $K\left(\rho \right)\lambda =-\frac{\partial J}{\partial u}$;
• Compute: ${▽}_{\rho }J=\frac{\partial J}{\partial \rho }+{\lambda }^T\left[\frac{\partial F\left(\rho \right)}{\partial \rho }-\frac{\partial K\left(\rho \right)}{\partial \rho }u\right]$.

3.3. Compliance Functional

The compliance can be defined as $C\left(\rho \right)=Fu\left(\rho \right)$, and its sensitivity can be determined as follows [97,98,99]:

$$\frac{\partial C}{\partial {\rho }_e}=F\frac{\partial u}{\partial {\rho }_e};$$

(43)

since $u\left(\rho \right)=K^{-1}\left(\rho \right)F$ and $\frac{\partial u}{\partial {\rho }_e}=-K^{-1}\frac{\partial K}{\partial {\rho }_e}{K}^{-1}F$, Equation (43) becomes:

$$\frac{\partial C}{\partial {\rho }_e}=-F{K}^{-1}\frac{\partial K}{\partial {\rho }_e}{K}^{-1}F=-u\frac{\partial K}{\partial \rho }u.$$

(44)

Considering that $K=A{\int }_{{\mathsf{\Omega }}_e}{B}_a^T\mathbb{C}{B}_{a}dx$, with A being the assembly operator, and that $\rho$ is constant in ${\mathsf{\Omega }}_e$, $\frac{\partial K}{\partial \rho }$ can be rewritten as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial K}{\partial \rho }& ={\int }_{{\mathsf{\Omega }}_e}{B}_a^T{\mathbb{C}}^{\prime }{B}_{a}dx={\int }_{{\mathsf{\Omega }}_e}{B}_a^{T}P{\rho }^{P-1}{\mathbb{C}}_0{B}_{a}dx={\int }_{{\mathsf{\Omega }}_e}\frac{P}{\rho }{B}_a^T\mathbb{C}{B}_{a}dx=\frac{P}{\rho }K.\hfill \end{array}\end{array}$$

(45)

Introducing Equation (45) into (44) leads to:

$$\begin{array}{c}\hfill \frac{\partial C}{\partial {\rho }_e}=-\frac{P}{\rho }{\rho }^P{u}_e^T{K}_0{u}_e=-P\frac{E_e}{\rho }\end{array}$$

(46)

where $u_e$ and $K_0$ are the elemental displacement vector and stiffness matrix of a fully solid element (i.e., $\rho =1.0$). The term ${\rho }^P{u}_e^T{K}_0{u}_e$ is the strain energy $\left(E_e\right)$, missing only the $\frac{1}{2}$ constant. However, because this constant is applied to all elements, it can be neglected, and set the term ${\rho }^P{u}_e^T{K}_0{u}_e$ equal to the strain energy ($E_e$) automatically calculated in ABAQUS^®.

In this particular case, the first three steps of the generic optimization process can be simplified and rewritten as follows:

• Find u such that: $K\left(\rho \right)u=F\left(\rho \right)$. This procedure can be performed using an ABAQUS^® FEA;
• Take $\lambda =-u$, since the problem is self-adjoint;
• Compute: $\frac{\partial C}{\partial {\rho }_e}=-P\frac{E_e}{\rho }$. Note that $E_e$ can be obtained from the ABAQUS^® FEA executed in the first step.

This information is included to allow an easier understanding of the code implementation and its correlation with the formal mathematical formulation.

3.4. Stress Functional

The maximum function is approximated by means of a modified p-norm function. This approximation is necessary to provide a derivable function that approximates the maximum function, which is non-differentiable. This implementation adopts the modified p-norm approximation proposed in [53]:

$${\sigma }^{PN}\left(\rho \right)={\left(\frac{1}{N_i}\sum _{\mathsf{\Omega }}{\left({\sigma }_a^{vM}\left(\rho \right)\right)}^q\right)}^{\frac{1}{q}}$$

(47)

where q is the exponential factor, $\mathsf{\Omega }$ is the set of stress evaluation points in the topology optimization problem, and ${\sigma }_a^{vM}$ is the value of the amplified von Mises stress in point a. Note that, although unconventional, this manuscript adopted q as the exponential factor of the p-norm approximation to improve the readability of the equation and avoid confusion with the design density symbol $\rho$ and the SIMP exponential factor P. Furthermore, it is relevant to highlight that, by adopting the use of this approximation, a global approach is selected, instead of local enforcement at each integration point. This option was preferred for providing better scaling with the mesh refinement, reducing the computational cost, and for its simplicity. However, the use of augmented Lagrangian methods could be seen as an alternative [57].

The derivative of the p-norm approximation with respect to the design density of an element can be obtained by the chain-rule, multiplying three intermediate terms. The first term is the derivative of the p-norm approximation with respect to the amplified von Mises stress:

$$\begin{array}{c}\hfill \frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}={\left(\frac{1}{N_i}\sum _{\alpha \in \mathsf{\Omega }}{\left({\sigma }_a^{vM}\left(\rho \right)\right)}^q\right)}^{\frac{1}{q}-1}\times \frac{1}{N_i}{\left({\sigma }_a^{vM}\left(\rho \right)\right)}^{q-1}.\end{array}$$

(48)

The second term is the derivative of ${\sigma }_a^{vM}$ with respect to the stress vector in point a. Since ${\sigma }_a^{vM}$ can be written in matrix form as ${\sigma }_a^{vM}={\left({\sigma }_{a}M{\sigma }_a\right)}^{\frac{1}{2}}$, its derivative becomes:

$$\begin{array}{c}\hfill \frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a\left(\rho \right)}={\left({\sigma }_a\left(\rho \right)M{\sigma }_a\left(\rho \right)\right)}^{-\frac{1}{2}}\times {\sigma }_a\left(\rho \right)M\frac{\partial {\sigma }_a\left(\rho \right)}{\partial \rho }.\end{array}$$

(49)

The third and last terms are the derivative of the stress vector ${\sigma }_a$ with respect to the design density. Considering Equations (9) and (10):

$$\begin{array}{c}\hfill \frac{\partial {\sigma }_a\left(\rho \right)}{\partial \rho }=\beta {\rho }^{\beta -1}{\mathbb{C}}_0{B}_{a}u\left(\rho \right)+{\rho }^{\beta }{\mathbb{C}}_0{B}_a\frac{\partial u\left(\rho \right)}{\partial \rho }.\end{array}$$

(50)

The term $\frac{\partial u\left(\rho \right)}{\partial \rho }$ is obtained from the state Equation (34):

$$\begin{array}{c}\hfill \frac{\partial K\left(\rho \right)}{\partial \rho }u\left(\rho \right)+K\left(\rho \right)\frac{\partial u\left(\rho \right)}{\partial \rho }=\frac{\partial F\left(\rho \right)}{\partial \rho }=0\iff \\ \hfill \iff \frac{\partial u\left(\rho \right)}{\partial \rho }=-K^{-1}\left(\rho \right)\frac{\partial K\left(\rho \right)}{\partial \rho }u\left(\rho \right)\end{array}$$

(51)

and, therefore, with $\beta =\frac{1}{2}$:

$$\begin{array}{c}\hfill \frac{\partial {\sigma }_a\left(\rho \right)}{\partial \rho }=\frac{1}{2}{\rho }^{-\frac{1}{2}}{\mathbb{C}}_0{B}_{a}u\left(\rho \right)-{\rho }^{\frac{1}{2}}{\mathbb{C}}_0{B}_a{K}^{-1}\left(\rho \right)\frac{\partial K\left(\rho \right)}{\partial \rho }u\left(\rho \right).\end{array}$$

(52)

Note that it was assumed that $\frac{\partial F\left(\rho \right)}{\partial \rho }=0$. This assumption is valid for static-load-driven problems, where the load is considered to be constant and independent of the material distribution. For displacement-driven problems, this assumption is not valid, as the forces resulting from the displacement applied will change depending on the material distribution.

With Equations (48) through (51) and applying the chain rule, it is possible to define $\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial \rho }$ as:

$$\begin{array}{c}\hfill \frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial \rho }=\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}\frac{\partial {\sigma }_a\left(\rho \right)}{\partial \rho }=\\ \hfill =\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}\left(\frac{1}{2}{\rho }^{-\frac{1}{2}}{\mathbb{C}}_0{B}_{a}u\left(\rho \right)\right.\\ \hfill \left.-{\rho }^{\frac{1}{2}}{\mathbb{C}}_0{B}_a{K}^{-1}\left(\rho \right)\frac{\partial K\left(\rho \right)}{\partial \rho }u\left(\rho \right)\right).\end{array}$$

(53)

Renaming these two terms as $\partial {\sigma }_{spf}^{PN}\left(\rho \right)$ (Equation (54)), referring to the component of the derivative that is dependent on the stress penalization factor, and $\partial {\sigma }_u^{PN}\left(\rho \right)$ (Equation (55)), referring to the component of the derivative that is dependent on the nodal displacement:

$$\begin{array}{c}\hfill \partial {\sigma }_{spf}^{PN}\left(\rho \right)=\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}\left(\frac{1}{2}{\rho }^{-\frac{1}{2}}{\mathbb{C}}_0{B}_{a}u\left(\rho \right)\right)\end{array}$$

(54)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \partial {\sigma }_u^{PN}\left(\rho \right)& =\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}\times \left(-{\rho }^{\frac{1}{2}}{\mathbb{C}}_0{B}_a{K}^{-1}\left(\rho \right)\frac{\partial K\left(\rho \right)}{\partial \rho }u\left(\rho \right)\right).\hfill \end{array}\end{array}$$

(55)

Notice that $\partial {\sigma }_{spf}^{PN}\left(\rho \right)$ can be easily determined, since the only information required is the stress vector $({\widehat{\sigma }}_a\left(\rho \right)={\mathbb{C}}_0{B}_{a}u\left(\rho \right))$ and the design densities. $\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}$ can be determined with a simple summation and $\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}$ through the product of two vectors and a matrix.

The term $\partial {\sigma }_u^{PN}\left(\rho \right)$, on the other hand, requires the explicit definition of the matrices $B_a$, $K^{-1}\left(\rho \right)$, and $\frac{\partial K\left(\rho \right)}{\partial \rho }$, which are dependent on the element formulation used in the numerical model. Furthermore, for topology optimization problems where the number of design variables is larger than the number of constraints, the most-efficient way to determine $\partial {\sigma }_u^{PN}\left(\rho \right)$ is using an adjoint model, defining the adjoint variable as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill K\left(\rho \right)\lambda & =\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}{\mathbb{C}}_0{B}_a\iff \hfill \\ \hfill \iff \lambda & =\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}{\mathbb{C}}_0{B}_a{K}^{-1}\left(\rho \right).\hfill \end{array}\end{array}$$

(56)

Therefore, the adjoint variable can be extracted from a finite element analysis, where the load applied on each node is equal to $\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}{\mathbb{C}}_0{B}_a$, and $\lambda$ is equal to the node displacement. Introducing the adjoint variable in Equation (54) then leads to:

$$\begin{array}{c}\hfill \partial {\sigma }_u^{PN}\left(\rho \right)=\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}\left(-{\rho }^{\frac{1}{2}}\lambda \frac{\partial K\left(\rho \right)}{\partial \rho }u\left(\rho \right)\right).\end{array}$$

(57)

Finally, since K can be defined as $K=A{B}_a^T\mathbb{C}\left(\rho \right)B_a$, since $B_a\lambda$ and $B_{a}u$ are equal to the deformation vectors of the adjoint and state models (${\xi }_{adj}$ and ${\xi }_s$, respectively), Equation (57) can be simplified to:

$$\begin{array}{c}\hfill \partial {\sigma }_u^{PN}\left(\rho \right)=\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}}{\partial {\sigma }_a}\left(-{\rho }^{\frac{1}{2}}{\xi }_{adj}^T\left(\rho \right){\mathbb{C}}^{\prime }\left(\rho \right){\xi }_s\left(\rho \right)\right).\end{array}$$

(58)

With this information, the first three steps of the generic optimization process can be rewritten as follows:

• Find u such that: $K\left(\rho \right)u=F\left(\rho \right)$. This procedure can be performed through an ABAQUS^® FEA;
• Find $\lambda$ such that: $K\left(\rho \right)\lambda =\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}\left(\rho \right)}{\partial {\sigma }_a}{\mathbb{C}}_0{B}_a$;
• Compute: $\partial {\sigma }_u^{PN}\left(\rho \right)=\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_a^{vM}}\frac{\partial {\sigma }_a^{vM}}{\partial {\sigma }_a}\left(-{\rho }^{\frac{1}{2}}{\xi }_{adj}^T\left(\rho \right){\mathbb{C}}^{\prime }\left(\rho \right){\xi }_s\left(\rho \right)\right)$.

As demonstrated in Equation (55) through (58), determining $\partial {\sigma }_u^{PN}\left(\rho \right)$ is significantly more complex and computationally expensive, due to the necessity of using the adjoint model. Some researchers [100] have proposed approximating this term as $\partial {\sigma }_u^{PN}\left(\rho \right)=0$, in an attempt to reduce the complexity and computational cost of stress-dependent topology optimization problems. While doing so does have evident advantages, the results obtained for each element of the benchmark problem described in Section 7.2 show that, on average, $|\partial {\sigma }_u^{PN}\left(\rho \right)|>|\partial {\sigma }_{spf}^{PN}\left(\rho \right)|$. As a result, assuming that $\partial {\sigma }_u^{PN}\left(\rho \right)=0$ would lead to an average error greater than $50\%$ in the calculation of the maximum stress sensitivity, at least for the particular case of the benchmark problem described in Section 7.2. The error in this approximation is aggravated by the fact that $\partial {\sigma }_u^{PN}\left(\rho \right)$ has a negative sign from Equation (53), and as a consequence, ignoring this term could potentially point the optimization algorithm in the opposite direction of the gradient. Therefore, the simplification proposed in [100] was not adopted in this work.

3.5. Volume Constraint

The value of the sensitivity of the volume constraint to changes in the design density of each element is equal to the value of its volume, as defined in Equation (59). Note that the equivalent sensitivity can be obtained for a mass constraint, replacing the volume of the element with its mass.

$$\frac{\partial V\left(\rho \right)}{\partial {\rho }_e}=V_e.$$

(59)

4. Mesh Dependency and Data Filtering

In order to obtain a mesh-independent solution and avoid the “checkerboard” instability, the raw sensitivities and/or design densities are processed with a blurring filter [6,101]. This article adopted the filter scheme used in [17], which is a simplification of the scheme proposed by Huang and Xie [102], described as follows:

$${\rho }_e=\sum _j\left(\frac{w\left(r_{ej}\right)}{{\sum }_{j}w\left(r_{ej}\right)}{\rho }_j\right)$$

(60)

where the value of $w\left(r_{ej}\right)$ is equal to the difference between the maximum range of the filter $\left(r_{max}\right)$ and the actual distance between the central element e and the j elements in its neighborhood $\left(r_{ej}\right)$, defined as follows:

$$w\left(r_{ej}\right)=max(0,r_{max}-r_{ej}).$$

(61)

This parameter can be interpreted as a pondered measurement of how close the two elements are. Note that the code implementation provided allows the user to select if the blurring filter should be applied to the sensitivities (excluding the volume constraint sensitivities, or equivalent), the design densities, or both.

It is relevant to note that, although not implemented in the methodology herein proposed, it is common to couple filtering techniques with projection methods [103] in order to obtain a better geometry definition in the final solution. Doing so may eliminate residual grey areas, which is particularly relevant in stress-constrained problems, where a smoother geometry may avoid the creation of stress concentration points.

5. Optimization Algorithms

5.1. Optimality Criteria

The OC method implemented in this article follows the approach proposed by Bendsøe [99] and implemented by Sigmund in [6]. It is usually applied to problems of compliance minimization with a volume constraint. According to this procedure, the design variables can be updated according to the following expression:

$$\rho =\left\{\begin{array}{cc}{\rho }_{lower},\hfill & \mathrm{if}\rho B_e^{\eta }\le {\rho }_{lower}\hfill \\ \rho B_e^{\eta },\hfill & \mathrm{if}{\rho }_{lower}\le \rho B_e^{\eta }\le {\rho }_{upper}\hfill \\ {\rho }_{upper},\hfill & \mathrm{if}{\rho }_{upper}\le \rho B_e^{\eta }\hfill \end{array}\right.$$

(62)

$${\rho }_{lower}=max({\rho }_{min},\rho -\delta \rho )$$

(63)

$${\rho }_{upper}=min(1.0,\rho +\delta \rho )$$

(64)

where ${\rho }_{lower}$ and ${\rho }_{upper}$ represent the design densities’ move limits imposed by the parameter $\delta \rho$. $\eta$ is a numerical damping coefficient set equal to $0.5$ in accordance with [6]. $B_e$ is a parameter obtained from the optimality condition, described as the ratio between the gradient of the objective and constraint functions divided by the Lagrangian multiplier ($\psi$), as follows:

$$B_e=\frac{-\frac{\partial C}{\partial \rho }}{\psi \frac{\partial V}{\partial \rho }}=\frac{-\frac{\partial C}{\partial \rho }}{\psi V_e}.$$

(65)

The Lagrangian multiplier $\psi$ can be found using a bisection method, updating its value until the volume of the solution meets the imposed constraint.

This method can also be simplified in order to obtain an equivalent version suitable for discrete design variables. This simplification corresponds to the approach used by Zuo and Xie in [17]. For the current iteration, the code determines the sensitivities of all elements. Then, a bisection method is used to determine a threshold sensitivity. All elements with a sensitivity greater than the threshold become solid (i.e., $\rho =1.0$), while the other elements are removed according to the “soft-kill” approach (i.e., $\rho ={\rho }_{min}$). Similar to the continuous version, the value of the sensitivity threshold is updated until the mass or volume of the structure meets the volume constraint imposed.

When using the OC method, especially the discrete version, it is usual to consider a fully solid design for the initial iteration and gradually reduce the volume constraint in each iteration until the target value is reached. This procedure is adopted to obtain a convergent solution and can be defined through the following equation:

$$V^{k+1}=max(V^{*},V^k(1.0-e_{vol}))$$

(66)

where k is the number of the current iteration and $e_{vol}$ is the ratio at which the volume constraint decreases in each iteration, until the target volume constraint ($V^{*}$) is reached.

Nonetheless, it is important to understand that this discrete version should be regarded as a heuristic, which was observed to work well when solving compliance-minimization problems. However, this procedure can be criticized from an optimization perspective, as there is no mathematical guarantee that the objective function will decrease during each iteration.

5.2. Method of Moving Asymptotes

The Method of Moving Asymptotes (MMA) is a nonlinear programming method that can be applied to structural optimization problems. Its functioning is based on an iterative process, generating a series of strictly convex approximating sub-problems. As an input to generate the sub-problem approximation, the MMA takes the value of the objective function, the value of the constraints, and the respective gradients as the inputs. Then, using dual-methods such as [104,105], the sub-problem is solved and the solution is taken as the value of the design variables in the next iteration of the main optimization problem [82].

The popularity of the application of the MMA to topology optimization problems is justified by several factors. From a computational point of view, the use of an approximation sub-problem and gradient information provides a means to reduce the number of structural problems solved, which are usually expensive. Furthermore, the versatility of the method and capability of handling multiple types of constraints make it suitable and applicable to a wide range of cases, from simple compliance-minimization problems with a single volume constraint, to cases where multiple constraints with different natures are applied. Its stability is also a positive characteristic, even for initial non-fully solid designs. Due to its extension and complexity, the interested reader is referred to [82] for a detailed explanation of this method.

The present code establishes a connection between ABAQUS^® and the implementation of the MMA created by Deetman [106] (“mmasub” and “subsolve” functions of the code provided), which is a Python version of the MATLAB code created by Svanberg [107].

5.3. Sequential Least-Squares Programming

Sequential Least-Squares Programming (SLSQP) is a Sequential Least-Squares Programming algorithm that utilizes the Quasi-Newton method to obtain an optimal solution. The strategy behind this algorithm consists of approximating the objective function with a quadratic equation and using the minimum of the approximation function as a possible solution in the next iteration. To improve the efficiency, the SLSQP can use Quasi-Newton methods to approximate the objective Hessian function.

The code provided uses the SLSQP implementation available in the SciPy module [108], which is a library of numerical routines for the Python programming language, which provides the fundamental building blocks for modeling and solving scientific problems.

5.4. Trust-Constr

Trust-constr is a trust region algorithm for constrained optimization, available in the SciPy module [108]. Both line search methods and trust region methods generate steps with the help of a quadratic model of the objective function. However, while line search methods use a step length along a given search direction, trust region methods define a region around the current iteration, within which they trust that the quadratic model is an accurate approximation [109].

This algorithm utilizes the trust region interior point method described in [110] to solve the inequality constraints imposed in the topology optimization problem statement. The strategy behind this interior point method [110] consists of solving inequality constraints by introducing slack variables. Then, the process is repeated for a sequence of equality-constrained barrier problems with progressively smaller values of the barrier parameter.

6. Python Implementation and Usage

The following subsections present a brief guide on how to use the Python code provided with the following Dataset [94] or repository: https://github.com/pnfernandes/Python-Code-for-Stress-Constrained-Topology-Optimization-in-ABAQUS. A brief overview of the most-important details of the main code classes and functions is also included in these subsections. The detailed information of each class and function can be found in their respective doc-string. The sequence of operations performed by this code implementation is described in Figure 1.

This code can be called in ABAQUS^® using the Run Script command in the File tab or by copying it into the command line. Note that the Run Script command loads the code faster, but automatically assumes that the user intends to solve the topology optimization in its totality. Copying it into the command line leads to a slower input, but allows a line-by-line execution, which may be useful for understanding the functioning of the code implementation.

The PEP8 [111] Python style guide was followed as closely as possible. However, the interested reader is warned that some exceptions exist, caused by compatibility constraints with the data structures and pre-defined variables existing in ABAQUS^®. An example of these exceptions includes the naming of constant material properties without all capital letters, since these names are already used by ABAQUS^® in a different data structure.

6.1. Code Usage

Upon executing the code, ABAQUS^® will create a series of prompt dialogue boxes. This allows the user to select the model database file (.cae), structural part, material, problem statement, optimization algorithm, and internal parameters to be considered during the topology optimization.

Besides being fully functional, the ABAQUS^® model does not require any specific formatting, except for the following considerations:

• Within the same part, the user can define specific elements to be optimized, dividing the part into two regions: one with editable elements and another with elements that should not be included in the optimization process (i.e., passive elements, from this point onward referred to as “frozen region” or “frozen elements”). To do so, the user should create a set named “editable_elements”, containing the elements to be edited. If this set does not exist, the whole part will be optimized. Note that the name “editable_elements” is case sensitive. This is performed in order to promote the optimization only of the targeted area and also to reduce the computational cost.
• The region to be optimized should not have a material section assigned to it, as the code will automatically assign the material properties and sections directly to the elements. Therefore, if the user defined a set with the editable elements, only the elements not included (frozen elements) should have a material section assigned to them.
• It is possible to allow the frozen elements to be considered during the filtering operations. To do so, the user should create a second element set with the case-sensitive name “neighbouring_region”.
• If multiple copies of the part to be optimized are included in the model assembly, only the first copy will be optimized.

6.2. Model Formatting, Job Submission, and Sensitivities (Lines 28–3452)

The ModelPreparation class modifies the ABAQUS^® model file (.cae) in order to allow the assignment of material properties and material sections to each individual element, as well as the property update during each iteration. Furthermore, this class requests the output of the strain energies, external work, and node strains if solving a stress-dependent problem. In ABAQUS^®, the strain energies of geometrically nonlinear problems are stored in the variables “SENER” and “PENER”, respectively, for the elastic and plastic strain components. Note that, according to the literature [14,17,98,102,112], the sensitivity of each element in a geometrically nonlinear compliance-minimization problem is given by the sum of the elastic strain $\left(E_e^e\right)$ and plastic strain components $\left(E_e^p\right)$, as described in Equation (67). However, for geometrically linear problems, the plastic strain component is zero, and ABAQUS^® stores the strain energies in the variable “ESEDEN”. The external work is only requested for a critical analysis of the result, as it should be equal to the objective function considering no energy dissipation.

$$\frac{\partial C}{\partial \rho }=E_e^e+E_e^p.$$

(67)

The submission of the FEA simulations for the state and adjoint models is managed by the classes ABAQUSFEA and AdjointModel, respectively. ABAQUSFEA is also responsible for determining the sensitivity of the compliance objective function $\left(\frac{\partial C}{\partial \rho }\right)$, which is stored as the class attribute “ae”. On the other hand, the AdjointModel class also applies the adjoint loads and determines the stress sensitivity. The stress sensitivity $\left(\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial \rho }\right)$, and its most-relevant intermediate terms ($\partial {\sigma }_{spf}^{PN}\left(\rho \right)$ and $\partial {\sigma }_u^{PN}\left(\rho \right)$) are stored as attributes of this class, along with the stress and strain vectors. In particular:

• “elmt_stress_sensitivity_continuous” stores the values of $\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial \rho }$ in a mesh-independent and mesh-independent form, making the comparison of the value of this derivative more suitable for the results obtained through finite differences.
• “elmt_stress_sensitivity_discrete” stores the values of $\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial \rho }$ in a mesh-dependent form, which is more suitable to be used as an input to the optimization algorithms, since they generally do not have access to the information of the mesh used in the FEA model. The difference between discrete and continuous element stress sensitivity is the multiplication by the determinant of the Jacobian matrix.
• “d_pnorm_spf” stores the values of $\partial {\sigma }_{spf}^{PN}\left(\rho \right)$.
• “d_pnorm_displacement” stores the values of $\partial {\sigma }_u^{PN}\left(\rho \right)$.
• “stress_vector” and “stress_vector_int” store the stress vectors determined at the element nodes and integration points.
• “deformation_vector” and “deformation_vector_int” store the strain vectors determined at the element nodes and integration points.

The sensitivity of the material constraint is determined by the material_constraint_sensitivity function, accounting for the possibility of having a non-uniform mesh with elements of different sizes and the use of either a volume or mass constraint.

6.3. Material and Stress Constraints (Lines 3453–3607)

The values of the material and maximum stress constraints are determined by the MaterialConstraint class and the stress_constraint_evaluation function, respectively. The maximum stress value is approximated with the modified p-norm function described in Equation (47), which is determined by the p_norm_approximation function.

6.4. Data Filtering (Lines 3608–3846)

The blurring filter is defined by the DataFilter class, which determines how each element is influenced by its neighborhood. The neighborhood of each element is defined by the elements that are fully within a maximum search radius defined by the user. The code Lines 3760–3762, intentionally left commented to reduce computational cost, generate element sets that allow an easy visual interpretation of the neighborhood of each element.

6.5. Optimization Algorithms: Optimality Criteria, Method of Moving Asymptotes, Sequential Least-Squares Quadratic Programming, and Trust-Constr (Lines 3847–5622)

The present code provides five optimization algorithms, described in Section 5, which can be used to solve the topology optimization problem.

The discrete and continuous versions of the OC algorithm, implemented in the functions oc_discrete and oc_continuous, are only suitable for compliance-minimization problems. The remaining optimization algorithms are suitable for all problem statements included.

The mma function is a decorator that links ABAQUS^® to the mmasub and subsolv functions implemented by Arjen Deetman [106]. Therefore, the main purpose of the mma function is the processing of the inputs necessary to apply the MMA [107] as a function of the problem statement selected by the user.

The functioning of the oc_discrete, oc_continuous, and mma functions is managed by the main program loop, described in Section 6.11. However, this is not the case for the SLSQP and Trust-constr algorithms available in this code, as they belong to the SciPy module [108]. To be fully functional, the SLSQP and Trust-constr algorithms require the additional freedom to decide when and how many times to call the objective, sensitivity, and constraint functions. Furthermore, the problem statement must be defined using unique data structures. For these reasons, class SciPyOptimizer is responsible for:

• Calling the SLSQP and Trust-constr methods from the SciPy module [108].
• Re-defining the problem statement with the correct data structure required by the algorithm selected.
• Providing the freedom for these algorithms to call the objective, sensitivity, and constraint functions as many times and in any given order necessary.
• Trying to record the data of each iteration as accurately as possible. However, it is to be expected that the data recorded may include information from intermediate evaluation points.

Finally, there are two details that should be considered before using the SLSQP and Trust-constr algorithms. One is that both have internal convergence criteria. To change them, the user may be required to edit the code provided with this article. The second detail is that the initial solution must be feasible; otherwise, the optimization process will stop. This is particularly relevant when defining the initial material distribution and/or when using a p-norm continuation approach (described in Section 6.11). The OC and MMA methods implemented do not have these restrictions for the initial guess.

6.6. Display Definition (Lines 5623–6207)

The SetDisplay class is capable of changing the color codes of the element sets. In compliance-minimization problems, this allows the user to plot a gray-scale representation of the material distribution. In stress-dependent problems, the user has the additional options of plotting the distribution of the element stress or element-amplified stress, both of which can be raised to the p-norm exponential factor. The element stress plotted is equal to the average of the stresses determined at the integration points of the element.

The graphic representation(s) plotted at each iteration are automatically saved to a .png file.

To allow an easier interpretation of the 3D material distributions, the SetDisplay class has a built-in method (function of a class) named hide_elements, whose input is a design variable threshold value. This method will hide all elements whose design variable is below the threshold input.

6.7. Data Recording (Lines 6208–6349)

At the end of each iteration, the save_data function will create a text file that records the values of all variables necessary to describe the current solution, as well as the internal parameters selected by the user.

This text file can also be used to restart the topology optimization process, starting from any iteration. However, as mentioned in Section 6.5, note that the algorithms SLSQP and Trust-constr may require the current iteration to be feasible with respect to the constraints imposed.

At the end of the optimization process, save_mdb will create a separate ABAQUS^® .cae file with the final solution obtained and a record of the objective function and material constraint values.

6.8. Element Formulation and Stiffness Matrix (Lines 6350–7283)

To solve stress-dependent problems, it is necessary to have access to information that depends on the element formulation (as shown in Section 3.4), such as $B_a$, K, $\mathbb{C}$, and the Jacobian matrix, as well as the definition of the element shape functions and their derivatives. The code provided includes the information required to use one 2D element with 4 nodes and one 3D element with 8 nodes. The 2D element is referred in ABAQUS^® by the code 2DQ4 and, more specifically, as CPS4 or CPE4 for the plane–stress and plane–strain cases, respectively. The 3D element is referred in ABAQUS^® by the code C3D8. The information of these elements is included in the ElementFormulation class, while the stiffness matrix $\mathbb{C}$ is created by the c_matrix_function function.

Although not used in the present article, the ElementFormulation class also includes the formulation of a shell element (referenced by the code S4 in ABAQUS^®). This information was included in order to promote the development of topology optimization methodologies suitable to this particular type of element.

6.9. Parameter Input Request, Domain Definition, and Variable Generation (Lines 7284–8854)

The ParameterInput class generates the prompt boxes that appear when running the code. The information introduced by the user is then used to create the global variables required for the code. A detailed list of these variables and their purpose can be found in Lines 8366–8461.

This information is then used by the EditableDomain class, to define the design space of the topology optimization problem, and by the VariableGenerator class, to create the lists and dictionaries that record the relevant data gathered during the topology optimization process. For a detailed description of these variables, please refer to the code Lines 8721–8827.

6.10. Auxiliary Functions (Lines 8855–9013)

Finally, there are four auxiliary functions to be briefly mentioned:

• average_ae outputs the average of the compliance sensitivity of each element over the last three iterations. This leads to an easier convergence of the OC algorithms, especially when considering discrete variables.
• The update_past_info function updates the variables that record the design variables and compliance sensitivity used in the previous two iterations.
• The evaluate_change function creates a ratio that describes how the objective function has changed over the last 10 iterations, which is used as convergence criteria in the implementation provided.
• The remove_files function will allow an automatic removal of the temporary files created by ABAQUS^® after each FEA.

6.11. Main Program (Lines 9014–9307)

The main program is divided in two phases. Between Lines 9014 and 9110, the code creates the classes required for the problem statement selected and prepares the ABAQUS^® model accordingly. The optimization process and the convergence criteria are defined in Lines 9111–9307.

The optimization process considers two loops. The first is applicable for stress-dependent problems, allowing a continuous increase of the p-norm exponential factor between $Qi$ and $QF$. This functionality, here referred to as the “p-norm continuation approach”, allows the stress constraint to be applied gradually and may be useful to improve the convergence of the algorithm. The continuation approach contrasts with the constant approach, which only considers one constant value for the p-normfactor. Note that for stress-independent problems (where $Qi$ is not used, but set to 1.0) and when using the constant approach, $Qi=QF$, leading to a single loop.

The second loop represents the convergence criteria. The code assumes that the algorithm has converged if the objective function has not changed more than 0.1 % over the last 10 iterations. It was selected as the default convergence criteria to allow an easy and relatively fast recreation of the results detailed in this article, as well as for being similar to the criteria adopted in [17]. The interested user is reminded that these criteria do not apply to the SLSQP and Trust-constr algorithms (as mentioned in Section 6.5) and is encouraged to create convergence criteria adapted to the topology optimization problem analyzed.

7. Benchmark Problems and Case Study

7.1. Cantilever Beam

The purpose of including this problem was to validate the functioning of the compliance objective function and volume constraint implemented in the code provided, before addressing more-complex stress-dependent problems, such as the L-bracket detailed in Section 7.2.

This problem was modeled using a regular mesh with an element size of 5.0 mm and the element type CPS4. The load F was equal to 100.0 N, applied to a single node, in the corner represented in Figure 2. The problem was solved using an implicit analysis and considering a material with a 70.000 MPa Young’s modulus, a 0.33 Poisson’s ratio, and a material density of $2.7\times {10}^{-9}$ t/mm³.

The ABAQUS^® models necessary to reproduce all the results described in this work are available in the Dataset [94].

7.2. L-Bracket

The L-bracket was chosen as a more-challenging benchmark problem [53,113,114], characterized by its initial geometry containing a stress-concentration point in the internal corner. Unconstrained compliance-minimization problems tend to ignore the stress concentration point, leading to an angular shape. On the other hand, stress-minimization or stress-constrained compliance-minimization problems tend to create rounded shapes that avoid the stress concentration. For these reasons, the L-bracket problem will be used in this research to validate and demonstrate the functioning of the code implemented when solving stress-dependent topology optimization problems. Furthermore, the strain and stress state generated by this L-bracket geometry is used in Section 8 to validate the correct implementation of the element formulations and maximum stress differentiation process.

The dimensions and boundary conditions considered are represented in Figure 3. Unless stated otherwise, this problem was modeled using a regular mesh with an element size of 3.0 mm and the element type CPS4. The load F was equal to 1500.0 N, with this value being distributed over the nodes included in the 12.0 × 12.0 corner represented in Figure 3. The elements within the 12.0 × 12.0 corner are not editable during the optimization process. The problem was solved using an implicit analysis and considering a material with a 70.000 MPa Young’s modulus, a 0.33 Poisson’s ratio, and a material density of $2.7\times {10}^{-9}$ ton/mm³.

7.3. Bonded Support

The problem described in Figure 4 was selected to provide a three-dimensional case study that is more representative of the complex geometries that can be found in engineering design problems.

The design space was modeled using an irregular mesh with an average element size of 1.25 mm and hexahedral C3D10, leading to an approximate total of 370.000 elements. A unit load was applied at a reference point attached to the connection region of the structure, while the surface of the adhesive joint was considered fixed. To minimize the computational cost, two planes of symmetry were considered along the XY and YZ planes, reducing the total number of elements. The problem was solved using an implicit analysis considering a material with: a 70.000 MPa Young’s modulus, a 0.33 Poisson’s ratio, and a material density of $2.7\times {10}^{-9}$ ton/mm³. The problem also considered a simplified representation of an adhesive bond, accounting only for the elastic properties (a 1.500 MPa Young’s modulus and a 0.3 Poisson’s ratio). Due to the complexity of the design, the interested reader is referred to the Mendeley Dataset [94] for the complete description of the geometry considered in this case study and shown in Figure 4.

The material elements representing the adhesive bond are identified as “frozen elements”. This excludes the possibility of editing the material distribution in this region, which would influence the resistance of the joint due to the reduction of the bonded area and excludes the contribution of these elements during the use of the blurring filter, which would be physically unreasonable as the adhesive and substrate are made of different materials. On the other hand, the first layer of elements of the support in contact with the adhesive was identified as a “neighboring region”, which does not allow changes in the material distribution of these elements, but allows the consideration of their influence in the determination of the gradients. However, unlike the previous case, the elements in this region were considered during the filtering process, as the material is the same as the one assigned to the editable region of the topology optimization process. The “frozen elements” and “neighboring region” are shown in Figure 4 in red and blue, respectively.

The ABAQUS^® models necessary to reproduce all the results described in this work are available in the Dataset [94].

8. Code Validation

The data necessary to solve compliance-minimization problems can be directly obtained from ABAQUS^®, which is a commercially certified software. However, stress-dependent problems require the explicit programming of the element formulation and of the derivation process of the maximum stress function. The purpose of the following subsections is to present a brief validation of the code implementation. Section 8.1 validates the element formulations implemented, comparing the strains and stresses determined by the code and by ABAQUS^® at each integration point. Section 8.2 compares the derivative of the maximum stress determined by the code with the derivative obtained through finite differences.

In this section, the L-bracket benchmark problem described in Section 7.2 is used. The reason behind this choice is that the L-bracket problem leads to a complex strain and stress state, allowing a more-challenging and -generic testing of elements over a wider variety of loading conditions than the Cantilever beam problem.

8.1. Validation of the Element Formulation

The analysis of the element formulation is based on the measurement of a pondered error, described as follows:

$$E_{rr}^i=\delta {\gamma }_i\frac{\gamma i}{\left|\gamma \right|}.$$

(68)

where $E_{rr}^i$ is the pondered error determined for the component i of a generic vector $\gamma$ and $\delta {\gamma }_i$ is the difference between the component determined by ABAQUS^® and its code counterpart. $E_{rr}^i$ is determined as the product of $\delta {\gamma }_i$ by the ratio between the component ${\gamma }_i$ and the magnitude of the vector $\gamma$.

This metric is preferred over the direct comparison of $\delta {\gamma }_i$, since the latter is highly influenced by floating point errors. One of the reasons for this influence is that, while ABAQUS^® operates internally with double-precision, the output received by the code is in a single-precision format. This difference reduces the number of decimal places considered in each value and could lead to disproportional artificial errors in the vector components that contribute the least to the magnitude of the vector or in vectors whose magnitude tends towards zero. The second reason is that $\delta {\gamma }_i$ tends to estimate large error values for vector components that have a significantly smaller magnitude when compared to the magnitude of the vector.

Table 2. Average pondered error between the elements implemented in the code and the ABAQUS^® output. The dimensions are represented in millimeters.

	CPS4		CPE4		S4		C3D8
	$E_{\mathrm{rr}}^i$ (%)	Int. Points $E_{\mathrm{rr}}^i$ < 1% (%)	$E_{\mathrm{rr}}^i$ (%)	Int. Points $E_{\mathrm{rr}}^i$ < 1% (%)	$E_{\mathrm{rr}}^i$ (%)	Int. Points $E_{\mathrm{rr}}^i$ < 1% (%)	$E_{\mathrm{rr}}^i$ (%)	Int. Points $E_{\mathrm{rr}}^i$ < 1% (%)
${\xi }_{11}$	0.05	99.65	0.90	63.39	0.49	97.29	0.44	96.37
${\xi }_{22}$	0.15	99.55	0.97	63.31	0.44	97.60	0.44	96.37
${\xi }_{33}$	-	-	-	-	-	-	0.44	96.37
${\xi }_{12}$	0.04	99.51	0.04	99.51	1.95	32.45	<0.01	100.0
${\xi }_{13}$	-	-	-	-	-	-	<0.01	100.0
${\xi }_{23}$	-	-	-	-	-	-	<0.01	100.0
${\sigma }_{11}$	0.00	99.88	1.85	31.74	0.61	95.89	1.45	44.81
${\sigma }_{22}$	0.15	99.65	2.00	63.41	0.61	97.29	1.45	96.36
${\sigma }_{33}$	-	-	-	-	-	-	1.45	44.82
${\sigma }_{12}$	0.01	99.55	0.01	99.55	0.77	77.44	<0.01	100.0
${\sigma }_{13}$	-	-	-	-	-	-	<0.01	100.0
${\sigma }_{23}$	-	-	-	-	-	-	<0.01	100.0

summarizes the average pondered error observed for each integration point of four variants of the L-bracket numerical model described in Section 7.2, each with a different element type. For the particular case of using a 3D element (element type C3D8 in ABAQUS^®), a thickness of 1.0 mm was considered for the L-bracket geometry. Since the largest average pondered error observed is equal to 2.0%, it can be concluded that the element formulations included in the code provided were implemented successfully for the case selected. Furthermore, for each element type, the table indicates the percentage of integration points with a pondered error less than 1.0%.

Although not included in the present work due to its extension, the interested reader is informed that the Dataset [94] associated with this research includes a detailed comparison of every integration point, considering both the pondered and regular error measurements.

8.2. Validation of the Maximum Stress Derivative

To validate the stress norm derivative of the code provided, its result was compared with the derivative obtained through finite differences. To do so, the design density of the elements located at points A, B, and C (shown in Figure 3) were changed individually, and the value of the maximum stress approximation was used to apply the finite differences. This procedure was repeated for the three elements, considering three different mesh sizes of 3.0 mm, 1.0 mm, and 0.5 mm, each increasing the order of magnitude of the total number of elements by 1, and different design density changes of 0.1, and 0.01. These density changes were chosen to conform with the maximum stiffness element difference allowable by ABAQUS^®. The results obtained are summarized in

Table 3. Analysis of the continuous maximum stress derivative determined by the Python code and the derivative obtained through finite differences.

Mesh Size	Number of Elements	p-Norm (MPa) @ $\mathit{\rho }$ = 1.0	Element Location	Element Label	$\mathbf{\Delta }\mathit{\rho }$	p-Norm (MPa) @ $\mathit{\rho }=1-\mathbf{\Delta }\mathit{\rho }$	Finite Differences	Continuous Derivative
0.5	102,400	314.27741	A	69,184	0.1	335.12679	−208.4939	−129.6819
					0.01	315.23590	−95.8493
			B	102,332	0.1	314.27740	0.0000	0.0001
					0.01	314.27741	0.0000
			C	63,841	0.1	314.28104	−0.0363	−0.0301
					0.01	314.27784	−0.0438
1	25,600	272.96255	A	17,269	0.1	272.87736	0.8519	0.4426
					0.01	272.95502	0.7534
			B	25,560	0.1	272.96249	0.0006	0.0004
					0.01	272.96255	0.0000
			C	23,901	0.1	272.97007	−0.0752	−0.0685
					0.01	272.96326	−0.0713
3	2889	221.08948	A	1953	0.1	232.36352	−112.7404	−65.2413
					0.01	221.48312	−39.3637
			B	2875	0.1	221.08923	0.0025	0.0029
					0.01	221.08941	0.0066
			C	2683	0.1	221.13180	−0.4232	−0.4054
					0.01	221.09332	−0.3842

.

It can be observed that the results obtained through finite differences are in good agreement with the derivative obtained by the Python code implemented, especially when the mesh size and the design density decrease.

Table 4. Average value of the stress derivative, and its components, observed in the L-bracket model at different design density values, considering three different mesh sizes.

Mesh Size	$\mathit{\rho }$	$\frac{\mathit{\partial }{\mathit{\sigma }}^{\mathit{PN}}\left(\mathit{\rho }\right)}{\mathit{\partial }\mathit{\rho }}$	$\mathit{\partial }{\mathit{\sigma }}_{\mathit{spf}}^{\mathit{PN}}\left(\mathit{\rho }\right)$	$\mathit{\partial }{\mathit{\sigma }}_{\mathit{u}}^{\mathit{PN}}\left(\mathit{\rho }\right)$	$\frac{|\mathit{\partial }{\mathit{\sigma }}_{\mathit{u}}^{\mathit{PN}}\left(\mathit{\rho }\right)|}{|\mathit{\partial }{\mathit{\sigma }}_{\mathit{spf}}^{\mathit{PN}}\left(\mathit{\rho }\right)|}$
0.5	1	−0.0088	0.0004	−0.0092	24.037
	0.5	−2304.3703	0.5426	−2304.9129	4248.295
	0.1	−90.9833	1.2132	−92.1965	75.996
1	1	−0.0305	0.0013	−0.0318	23.864
	0.5	−7949.4335	1.8849	−7951.3184	4218.445
	0.1	−313.8380	4.2147	−318.0527	75.462
3	1	−0.2140	0.0096	−0.2236	23.289
	0.5	−3.4697	0.1086	−3.5784	32.936
	0.1	−2206.1706	30.3676	−2236.5382	73.649

indicates the average values of $\frac{\partial {\sigma }^{PN}\left(\rho \right)}{\partial \rho }$ and its components $\partial {\sigma }_{spf}^{PN}\left(\rho \right)$ and $\partial {\sigma }_u^{PN}\left(\rho \right)$ observed in the L-bracket model as a function of the element mesh size and design density. These data indicated that the approximation proposed in [100] (described at the end of Section 3.4) should not be adopted for the present problem.

9. Topology Optimization Results

Figure 5 summarizes the results obtained for the compliance minimization of the Cantilever beam. Figure 6 summarizes the results obtained for the topology optimization of the L-bracket considering the different problem statements, plotting the objective function at each iteration and illustrating one of the final geometries obtained. Finally, Figure 7 shows the topology-optimized geometry of the bonded support. The color code adopted here displays elements with $\rho =0$ in white, intermediate densities in grey, and elements with $\rho =1.0$ in black. These results were not adimensionalized.

The results obtained for each problem statement and any particular choice of parameters are discussed in the following subsections. In common, these examples share the volume constraint $V^{*}=0.5$. In the particular case of the solution “OC—Discrete and decreasing”, the initial density was set to 1.0 and gradually reduced by a factor of $e_{vol}=0.05$, while the remaining cases considered an initial density equal to 0.5. The maximum move limit was equal to 0.2 for the OC or MMA algorithms and 1.0 for the SciPy algorithms SLSQP or Trust-constr, corresponding to their pre-default values. The blurring filter was applied to both the sensitivity and design densities with a maximum search radius of 8.0 mm for the L-bracket, 12.0 mm for the Cantilever beam, and 3.5 mm for the bonded support, except for the discrete version of the OC, which only applied the blurring filter to the sensitivity. The influence of the frozen elements was also considered during the application of the blurring filter. For stress-dependent problem statements, a p-norm exponential factor $Qi=8.0$ was used.

Overall, the Trust-constr and SLSQP algorithms tended to require a larger number of iterations to reach convergence. The reason for this difference was justified by their internal convergence criteria, as described in Section 6.11.

In the results shown, the existence of few elements with intermediate design densities is, generally, to be expected for two reasons. First, applying a blurring filter will cause a gradual transition between the solid and void regions, avoiding a full “black-and-white” solution. Second, as described in Section 6.11, the convergence criterion depends only on the values of the objective function and does not impose a strict restriction on the final solution being constituted only by solid elements. These two factors justify the apparent better performance of the discrete version of the OC algorithm, since it only uses solid elements and does not have a blurring filter applied to the element densities. However, it is most likely that this optimization method converged to a better local minimum.

It is also relevant to note that the mesh size had an influence on the use of elements with intermediate design densities. To completely avoid the existence of residual grey areas, it is recommended to implement projection methods [103], to define the final contour of the geometry, or to apply a final threshold filter to the design densities of each element. However, despite the potential influence of internal parameters, such as the filtering approach, mesh refinement, and the choice of the optimization algorithm and its hyperparameters, the results presented in both case studies are in line with the information reported in previous publications [53].

Regarding the computational efficiency of the implementation proposed here, it was seen that it requires an average of 50 s per iteration when executed with an Intel^® Core^TM i7-8750H CPU @ 2.20 GHz. However, the use of external packages, to apply the SLSQP and Trust-constr algorithms available in Python, led to an increase of this value to an average of 380 s per iteration. These values include the whole process required to perform an iteration, including: the update of the material properties, the preparation of the ABAQUS^® input file, model execution, data extraction from the FEA solution, the calculation of gradients, the application of the optimization process to update the design densities, and the plot of the new material distribution.

Finally, the Cantilever beam and L-bracket problems solved in these sections could also be replicated in 3D space. This is possible as the code provided includes the implementation of a 3D element (C3D8 in ABAQUS^®) and as the problem solving process does not require any particular adaptation when being converted to 3D space.

The interested reader is informed that the Dataset [94] contains all the solutions obtained by each optimization algorithm, at each iteration. Due to its large extension, the information presented in this section is limited to the most-relevant solutions.

9.1. Cantilever Beam: Compliance Minimization Results

The objective functions obtained by each of the different optimization algorithms were in good agreement, converging towards an approximately equal value of compliance. However, a different trend can be observed for the “OC—discrete and decreasing” case. In this particular case, the initial solution considers a fully solid design that must gradually reduce its mass to meet the volume constraint, causing a consequent reduction of the structural stiffness and an increase of the compliance value.

This observation is also in line with the volume constraint graphic, which shows a gradual decrease until the value 0.5 is reached. On the other hand, the remaining algorithms maintained the volume constraint within a feasible domain.

The final Cantilever geometry obtained was similar across all cases considered and illustrated, in Figure 5, with the solution found by the MMA. In this particular case, the elements with intermediate densities may also be the result of using a coarser mesh.

9.2. L-Bracket: Compliance Minimization Results

The solutions obtained for the compliance-minimization problem were similar across all the algorithms used, differing only in the use of intermediate design densities in some elements. The solution obtained by the continuous version of the OC is shown in Figure 6A as an example of this geometry.

The volume constraint was respected across all iterations. This behavior was also observed for the problem statements presented in Section 9.3 and Section 9.4. For this reason, the graphic representation of the constraints was not included in these sections.

9.3. L-Bracket: Stress-Constrained Compliance Minimization Results

The stress-constrained compliance-minimization problem statement considered a maximum allowable stress ${\sigma }^{*}=200.0$ MPa. Imposing this constraint caused the MMA to converge towards a curved geometry, shown in Figure 6B, avoiding the creation of a stress concentration point in the inner corner of the L-bracket. This transition is in agreement with the results presented in [53].

It is important to note that the transition from an angular geometry to a curved geometry (from case A to B of Figure 6) is dependent on several factors and not guaranteed to occur. First, the use of the modified p-norm approximation (Equation (47) proposed in [53]) may lead to an underestimate of the maximum stress installed. This fact may cause the optimization algorithm to overlook the influence of the stress concentration point and justifies the selection of ${\sigma }^{*}=200.0$ MPa, a relatively low maximum value for the stress constraint. Second, the algorithm selected may adopt a strategy that is more or less prone to explore and accept solutions that do not respect the imposed constraints. The solution obtained by the SLSQP algorithm had a slightly curved shape closer to the region where the load was applied, but did not avoid the internal angular geometry. The solution obtained by the Trust-constr region was equal to the one obtained in the compliance-minimization problem statement, as the implementation available in SciPy only determined the maximum stress sensitivity for the first iteration. Although this strategy reduced the computational cost, it did not allow a rigorous evaluation of the influence of each element in the maximum stress.

Finally, it is also important to note that the p-norm exponential factor had a significant influence on the result obtained, especially in the use of elements with intermediate densities.

9.4. L-Bracket: Stress Minimization Results

Defining a stress-minimization problem statement may also lead to a transition from an angular to a curved geometry. This is particularly evident in the solution obtained by the SLSQP algorithm, shown in Figure 6C. However, this transition did not occur with the method available in ABAQUS^®. As a result, the SLSQP and Trust-constr algorithms converged to solutions with maximum stresses of 315 MPa and 330 MPa, respectively, representing a 35% reduction when compared to the 518 MPa installed in the solution obtained with ABAQUS^®.

The data obtained with the SLSQP indicated an unstable behavior with large variability in the objective function. To improve readability, only the cumulative best solution obtained with this method was included in the plots referring to this problem statement. This variability indicates that the SLSQP was more likely to perform larger density re-distributions, which may be a result of the strategy adopted by the SciPy implementation of this algorithm, and a result of the influence of their internal parameters, which largely affect the solutions obtained.

Similar to the previous problem statement, the p-norm exponential factor had a significant influence on the result obtained, especially in the use of elements with intermediate densities. Also, although the MMA algorithm was removing material from the stress concentration area, further research should be performed to find appropriate tolerance parameters.

9.5. Bonded Support: Compliance Minimization Results

The application of the compliance-minimization problem statement to the bonded support allows the removal of material in the central region, leaving material in the connection between the adhesive joint and the elements sustaining the external loads and in the lower surface of the design domain. To improve the readability and comprehension of the geometry obtained, Figure 7 only displays the elements where $\rho =1.0$. The results herein shown also included the deformations installed in the topology-optimized design and the complete geometry obtained by mirroring the initial model along the XY and YZ planes.

10. Conclusions

This paper addressed the difficulties of applying novel topology optimization methods in the industry, where the black-box system common to most commercial software limits the access to relevant information. Therefore, a new open methodology for solving non-self-adjoint topology optimization problems in ABAQUS^® was presented.

The functioning of this method was validated through two procedures: first, by comparing the results of the sensitivity analysis with finite differences, which indicated a good correlation between both; then, by solving two benchmark problems with the following algorithms: Optimality Criteria, Method of Moving Asymptotes, Sequential Least-Squares Quadratic Programming, and Trust-constr algorithm. The results obtained, convergence, and influence of the problem statement chosen were in agreement with other state-of-the-art results. In the stress-minimization case study, having access to optimization algorithms unavailable in commercial software allowed the methodology herein presented a maximum stress reduction of 35% when compared to the solution obtained in ABAQUS^®.

An implementation of this method is provided in Python, working in conjunction with ABAQUS^®, in order to allow a direct application to new problems. The code implementation, as well as the ABAQUS^® models necessary to recreate the work reported in this paper are available as Supplementary Material and can also be downloaded in digital format from the following Dataset [94] or repository: https://github.com/pnfernandes/Python-Code-for-Stress-Constrained-Topology-Optimization-in-ABAQUS.

Therefore, this article promotes the advance of the state-of-the-art by providing the tools required to solve stress-dependent geometrically complex problems, as well as a modular path that opens the possibility of implementing novel methods in an industrial environment.