## 3 Discrete Formulation of Topology Optimization Problem Statements and Sensitivities
### 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}(\rho)=\left(\frac{1}{N_{i}}\sum_{\Omega}\Bigl{(}\sigma_{a}^{vM}( \rho)\Bigr{)}^{q}\right)^{\frac{1}{q}} \tag{47}\]  
where \(q\) is the exponential factor, \(\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:  
\[\frac{\partial\sigma^{PN}(\rho)}{\partial\sigma_{a}^{vM}}=\left(\frac{1}{N_{i }}\sum_{a\in\Omega}\Bigl{(}\sigma_{a}^{vM}(\rho)\Bigr{)}^{q}\right)^{\frac{1} {q}-1}\times\frac{1}{N_{i}}\Bigl{(}\sigma_{a}^{vM}(\rho)\Bigr{)}^{q-1}. \tag{48}\]The second term is the derivative of \(\sigma_{a}^{pM}\) with respect to the stress vector in point \(a\). Since \(\sigma_{a}^{pM}\) can be written in matrix form as \(\sigma_{a}^{pM}=(\sigma_{a}M\sigma_{a})^{\frac{1}{2}}\), its derivative becomes:  
\[\frac{\partial\sigma_{a}^{pM}(\rho)}{\partial\sigma_{a}(\rho)}=(\sigma_{a}( \rho)M\sigma_{a}(\rho))^{-\frac{1}{2}}\times\sigma_{a}(\rho)M\frac{\partial \sigma_{a}(\rho)}{\partial\rho}. \tag{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):  
\[\frac{\partial\sigma_{a}(\rho)}{\partial\rho}=\beta\rho^{\beta-1}\mathbb{C}_{ 0}B_{a}u(\rho)+\rho^{\beta}\mathbb{C}_{0}B_{a}\frac{\partial u(\rho)}{ \partial\rho}. \tag{50}\]  
The term \(\frac{\partial u(\rho)}{\partial\rho}\) is obtained from the state Equation (34):  
\[ \begin{split}\frac{\partial K(\rho)}{\partial\rho}u(\rho)& +K(\rho)\frac{\partial u(\rho)}{\partial\rho}=\frac{\partial F(\rho)}{ \partial\rho}=0\Leftrightarrow\\ &\Leftrightarrow\frac{\partial u(\rho)}{\partial\rho}=-K^{-1}( \rho)\frac{\partial K(\rho)}{\partial\rho}u(\rho)\end{split}  \tag{51}\]  
and, therefore, with \(\beta=\frac{1}{2}\):  
\[\frac{\partial\sigma_{a}(\rho)}{\partial\rho}=\frac{1}{2}\rho^{-\frac{1}{2}} \mathbb{C}_{0}B_{a}u(\rho)-\rho^{\frac{1}{2}}\mathbb{C}_{0}B_{a}K^{-1}(\rho) \frac{\partial K(\rho)}{\partial\rho}u(\rho). \tag{52}\]  
Note that it was assumed that \(\frac{\partial F(\rho)}{\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}(\rho)}{\partial\rho}\) as:  
\[