## 3 3D Orthotropic Stiffness Formulation
### Viscous Regularization  
Even when only the energy regularization described in the previous section is used, some issues may arise: (i) Loss of positive definiteness of the tangent stiffness; (ii) some damage localization; and (iii) some numerical convergence issues. To overcome these issues, a simple/efficient procedure is to implement viscous regularization.  
Damage localizes in a narrow band when the tangent stiffness tensor is not positive definite, and the numerical solution depends on the numerical discretization (decreasing the element size in the localized zone decreases the computed energy dissipated). Therefore, the structural response is not objective because it does not converge to a single solution with mesh refinement. Standard Newton-Rapshon approaches also have convergence issues when the tangent stiffness is not positively definite [52].  
As mentioned above, the viscous regularization approach can help to solve some of these convergence problems. This algorithm is built on an artificial Duvaut-Lions viscosity model [53]. More details on the theoretical background of this method can be found in the initial published work of [54]. This algorithm improves convergence by making the tangent stiffness matrix of the softening material positive definite for small enough time increments. In this regularization procedure, the incremental Equation (38) provides a definition for the viscous damage variable.  
\[\delta d_{v}=\frac{1}{\eta}(d-d_{v})\times\delta_{t} \tag{38}\]  
where \(d\) is the damage variable determined as previously explained, \(d_{v}\) is the regularized viscous damage variable, and \(\eta\) is the viscosity stabilization factor. The viscous and non-viscous stresses are written as (39) and (40), respectively.  
\[\{\sigma_{v}\}=[C_{d}^{v}]\{\varepsilon\} \tag{39}\]  
\[\{\sigma\}=[C_{d}]\{\varepsilon\} \tag{40}\]  
The damaged elastic matrix, \(C_{d}^{v}\), is computed using the viscous damage variables for each failure mode. The rate of convergence of the model in the softening regime is typically improved by using viscous regularization with a small value for the viscosity parameter (small compared to the characteristic time increment), without degrading the quality of the numerical solution. The fundamental notion is that the viscous system's solution relaxes to that of the inviscid case as \(\eta\to 0\), where \(t\) represents time. For the time integration of the internal variables, an algorithm must be implemented. The internal variables can be updated using a backward-Euler approach [55], resulting in the viscous damage variable (41), by merging it with the preceding expression (38).  
\[d_{v}^{t+\Delta t}=d_{v}^{t}+\Delta d_{v}^{t}\Rightarrow d_{v}^{t+\Delta t}=d_ {v}^{t}+\frac{1}{\eta}\left(d^{t+\Delta t}-d_{v}^{t+\Delta t}\right)\times \Delta_{t} \tag{41}\]Now, it is possible to calculate the viscous damage while keeping the implicit time integration (42), utilizing data from the past viscous damage and the current non-viscous damage.  
\[d_{v}^{t+\Delta t}=\frac{\eta}{\eta+\Delta t}d_{v}^{t}+\frac{\Delta t}{\eta+ \Delta t}d^{t+\Delta t} \tag{42}\]  
It is important to guarantee that the level of viscosity is not that high, which is performed indirectly by comparing the level of elastic damaged energy (43) and the viscous energy dissipation (44), in which the first must always override the second. If the time step \(\Delta t\) if very small, it is possible that the first term of the sum in (44) dominates the second, providing an increase in viscous energy dissipation, which might compromise the accuracy of the results.  
\[E_{d}^{t+\Delta t}=E_{d}^{t}+\frac{1}{2}\Big{[}\sum_{i=1}^{n}\Big{(}\sigma_{v, i}^{t+\Delta t}+\sigma_{v,j}^{t}\Big{)}\Delta\epsilon_{i}\Big{]} \tag{43}\]  
\[E_{v}^{t+\Delta t}=E_{v}^{t}+\frac{1}{2}\Big{[}\sum_{i=1}^{n}\Big{(}\sigma_{v, j}^{t+\Delta t}+\sigma_{v,j}^{t}-\sigma_{i}^{t+\Delta t}-\sigma_{i}^{t} \Big{)}\Delta\epsilon_{i}\Big{]} \tag{44}\]