## 2 Continuous Formulation of Topology Optimization Problem Statements and Sensitivities
### 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:  
\[J(\rho,u(\rho))=\left(\int_{\Omega}\left(\sigma_{a}^{VM}(\rho,u(\rho))\right) ^{q}dx\right)^{\frac{1}{q}} \tag{24}\]  
where \(q\) is the exponential factor, and the von Mises amplified stress norm is \(\sigma_{a}^{VM}=(\sigma_{a}M\sigma_{a})^{\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}l(\rho,u(\rho))(\tilde{\rho})\) is defined as:  
\[\begin{split}& D_{\rho}l(\rho,u(\rho))(\tilde{\rho})=\\ &\quad=\left(\int_{\Omega}\left(\sigma_{a}^{VM}(\rho,u(\rho)) \right)^{q}dx\right)^{\frac{1}{q}-1}\times q\int_{\Omega}\sigma_{a}^{VM}(\rho, u(\rho))^{q-1}D_{\sigma_{a}}\sigma_{a}^{VM}(D_{\rho}\sigma_{a}(\tilde{\rho}))\,dx \end{split} \tag{25}\]  
where:  
\[D_{\rho}\sigma_{a}(\tilde{\rho})=\beta\rho^{\beta-1}\mathbb{C}_{0}\bigtriangledown ^{s}u\tilde{\rho}. \tag{26}\]  
The derivative in the direction \(w\), which allows us to define the adjoint problem, is:  
\[\begin{split}& D_{u}l(\rho,u(\rho))(w)=\\ &\quad=\left(\int_{\Omega}\left(\sigma_{a}^{VM}(\rho,u(\rho)) \right)^{q}dx\right)^{\frac{1}{q}-1}\times\int_{\Omega}\sigma_{a}^{VM}(\rho,u( \rho))^{q-1}D_{\sigma_{a}}\sigma_{a}^{VM}(D_{u}\sigma_{a}(w))\,dx\end{split} \tag{27}\]  
where:  
\[D_{\sigma_{a}}\sigma_{a}^{VM}(D_{u}\sigma_{a}(w))=(\sigma_{a}M\sigma_{a})^{- \frac{1}{2}}\sigma_{a}MD_{u}\sigma_{a}(w) \tag{28}\]with:  
\[D_{u}\sigma_{a}(w)=\rho^{\beta}\mathbb{C}_{0}\bigtriangledown^{s}w. \tag{29}\]  
Finally, before determining the value of \(DJ(\rho,u(\rho))\tilde{\rho}\), it is required to find \(D_{\rho}a(\rho,u(\rho),\lambda)\) as follows:  
\[D_{\rho}a(\rho,u(\rho),\lambda)\tilde{\rho}=\int_{\Omega}\bigtriangledown^{s} \lambda\mathbb{C}^{\prime}(\rho)\bigtriangledown^{s}u\tilde{\rho}. \tag{30}\]  
As stated previously, note that \(\bigtriangledown^{s}\lambda\) and \(\bigtriangledown^{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\({}^{\beta}\) FEA. Also, note that, for the stress functional, \(D_{\rho}l(\lambda)=0\). Finally, the first three optimization steps can then be rewritten as:  
* Find \(u\), the solution of \(a(\rho,u(\rho)v)=l(\rho,v)\);
* Using Equation (27), find \(\lambda\) such that \(a(\rho,\lambda,w)=-D_{u}l(\rho,u(\rho))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(\rho))\tilde{\rho}=[D_{\rho}J(\rho,u(\rho))\ +\ D_{\rho}a(\rho,u(\rho),\lambda)] \tilde{\rho}\).