## 3 3D Orthotropic Stiffness Formulation  
The orthotropic formulation is described using ABAQUS standard notation for both stress and strain tensors [36]. The compliance and stiffness matrices are used to calculate the elastic constitutive relation for stress and strain (1). \[\varepsilon_{i}=H_{ij}\sigma_{j}\quad\left[H_{d}\right]^{-1}=\left[C_{d}\right] \quad\sigma_{i}=C_{ij}\varepsilon_{j} \tag{1}\]  
with:  
\[ \begin{array}{l}\sigma_{1}=\sigma_{11};\sigma_{2}=\sigma_{22};\sigma_{3}= \sigma_{33};\sigma_{4}=\tau_{12};\sigma_{5}=\tau_{13};\sigma_{6}=\tau_{23}\\ \varepsilon_{1}=\varepsilon_{11};\varepsilon_{2}=\varepsilon_{22};\varepsilon _{3}=\varepsilon_{33};\varepsilon_{4}=\tau_{12};\varepsilon_{5}=\tau_{13}; \varepsilon_{6}=\tau_{23}\end{array}  \tag{2}\]  
To present a 3D orthotropic stiffness that is both coherent and thermodynamically acceptable, it is necessary to assemble the damage compliance matrix \(H_{d}\) exclusively using the diagonal term (8) to account for the damage behavior [29; 37]. It is then possible to assemble the damaged stiffness matrix \(C_{d}\), by inverting the previous compliance matrix (10). It is important to state the relevance of \(\Gamma_{d}\) (13), because this term ensures that no transversal stress occurs when subjected to uniform longitudinal stress. Just as in the previous 2D Hashin formulation [15], the tension and compression damage \((d_{t},d_{c})\) are activated using the corresponding effective stress space (3). For this work, damage variables are associated with the fiber \(d_{f}\), the matrix \(d_{m}\), and the interlaminar \(d_{i}\). \[\hat{\sigma_{i}}=\frac{\sigma_{i}}{(1-d_{i})} \tag{3}\]  
where:  
\[d_{1}=d_{f};d_{2}=d_{m};d_{3}=d_{i};d_{4}=d_{s12};d_{5}=d_{s13};d_{6}=d_{s23}\]  
\[d_{f}=\left\{ \begin{matrix}d_{ft}&if&\partial_{1}\geq 0\\ d_{fc}&if&\partial_{1}<0\end{matrix} \right. \tag{4}\]  
\[d_{m}=\left\{ \begin{matrix}d_{mt}&if&\partial_{2}+\partial_{3}\geq 0\\ d_{mc}&if&\partial_{2}+\partial_{3}<0\end{matrix} \right. \tag{5}\]  
\[d_{i}=\left\{ \begin{matrix}d_{it}&if&\partial_{3}\geq 0\\ d_{ic}&if&\partial_{3}<0\end{matrix} \right. \tag{6}\]  
\[\left\{ \begin{matrix}d_{s12}=1-\Big{(}1-d_{ft}\Big{)}\Big{(}1-d_{fc}\Big{)}(1- d_{mt})(1-d_{mc})\\ d_{s13}=1-\Big{(}1-d_{ft}\Big{)}\Big{(}1-d_{fc}\Big{)}(1-d_{it})(1-d_{ic})\\ d_{s23}=1-(1-d_{mt})(1-d_{mc})(1-d_{it})(1-d_{ic})\end{matrix} \right. \tag{7}\]  
The main advantage is that, by design, this matrix verifies thermodynamic admissibility principles when using Gibbs free energy (14) [38; 39]. \[\left[H_{d}\right]=\left[ \begin{matrix}\frac{H_{11}}{1-d_{f}}&H_{12}&H_{13}&0& 0&0\\ H_{21}&\frac{H_{22}}{1-d_{m}}&H_{23}&0&0&0\\ H_{31}&H_{32}&\frac{H_{33}}{1-d_{i}}&0&0&0\\ 0&0&0&\frac{H_{44}}{1-d_{s12}}&0\\ 0&0&0&0&0&\frac{H_{66}}{1-d_{s13}}\end{matrix} \right. \tag{8}\]  
\[\left\{\begin{matrix}H_{11}=\frac{1}{E_{1}}\\ H_{22}=\frac{1}  
\[[C_{d}]=\Gamma_{d} \begin{bmatrix}C_{11}&C_{12}&C_{13}&0&0&0\\ C_{21}&C_{22}&C_{23}&0&0&0\\ C_{31}&C_{32}&C_{33}&0&0&0\\ 0&0&0&C_{44}&0&0\\ 0&0&0&0&C_{55}&0\\ 0&0&0&0&0&C_{66}\end{bmatrix}  \tag{10}\]  
with:  
\[\left\{ \begin{array}{l}C_{11}=\Big{(}1-d_{f}\Big{)}E_{1}[1-(1-d_{m})(1-d_{i} )v_{23}v_{32}]\\ C_{22}=(1-d_{m})E_{2}\Big{[}1-\Big{(}1-d_{f}\Big{)}(1-d_{i})v_{13}v_{31}\Big{]} \\ C_{33}=(1-d_{i})E_{3}\Big{[}1-\Big{(}1-d_{f}\Big{)}(1-d_{m})v_{12}v_{21}\Big{]} \\ C_{12}=C_{21}=\Big{(}1-d_{f}\Big{)}(1-d_{m})E_{1}[v_{21}+(1-d_{i})v_{31}v_{23}] \\ C_{13}=C_{31}=(1-d_{i})(1-d_{m})E_{1}[v_{31}+(1-d_{m})v_{21}v_{32}]\\ C_{23}=C_{32}=(1-d_{m})(1-d_{i})E_{2}\Big{[}v_{32}+\Big{(}1-d_{f}\Big{)}v_{31}v_ {12}\Big{]}\end{array} \right. \tag{11}\]  
\[\left\{ \begin{array}{l}C_{44}=(1-d_{s12})G_{12}/\Gamma_{d}\\ C_{55}=(1-d_{s13})G_{13}/\Gamma_{d}\\ C_{66}=(1-d_{s23})G_{23}/\Gamma_{d}\end{array} \right. \tag{12}\]  
\[\left\{\begin{array}{l}\Gamma_{d}=\Big{(} \begin{array}{l}1-\Gamma_{fm}v_{1 2}v_{21}-\Gamma_{mi}v_{23}v_{32}-\Gamma_{fi}v_{13}v_{31}-2\times\Gamma_{fmi}v_ {12}v_{23}v_{31}\end{array}