## 2 Cohesive Zone Model  
Cohesive elements are widely employed in the simulation of delamination propagation in composite materials. Cohesive elements model the interaction between adjacent layers by defining the traction-separation relationship along the delaminated interface. These elements capture the nonlinear behavior of the interface, allowing for accurate simulation of crack initiation and propagation. The behavior of cohesive elements is governed by the cohesive law [30], which defines the traction-separation relationship along the interface or crack. The cohesive law describes how the traction (force per unit area) between the surfaces changes as the separation (relative displacement) between them increases. Different models can be found in literature, but the simplest and most widely used employs a bilinear traction-separation law, as shown in Figure 1.  
At low levels of separation (\(\delta\)), the cohesive traction rapidly grows following the penalty stiffness (K), which represents the strong bonding between the laminae. As the separation increases, the stress reaches the interface strength of the material (\(\tau_{c}\)) at the damage initiation displacement (\(\delta^{0}\)) and the cohesive traction starts to decrease, indicating the progressive degradation of the interface and the opening of the crack. When the traction becomes zero, at the final displacement (\(\delta^{f}\)), the element completely fails to represent the  
Figure 1: Constitutive response for traction-separation cohesive elements. propagation of the crack. The degradation of the interface stress is driven by the damage variable, \(d\), which can be calculated using Equation (1):  
\[d=\frac{\delta^{f}\big{(}\delta-\delta^{0}\big{)}}{\delta\big{(}\delta^{f}- \delta^{0}\big{)}} \tag{1}\]  
The integral of the traction-separation curve represents the fracture toughness of the material (\(G_{\text{C}}\)). In a delamination propagation analysis, the length of the process zone (LPZ) [26], \(l_{pz}\), is defined as the distance between the crack tip and the location where the maximum cohesive stress occurs. The LPZ is a material property that depends on the response of the different material to the presence of a crack. The value of \(l_{pz}\) can be estimated for mode I loading conditions as follow:  
\[l_{pz}=\gamma\frac{EG_{\text{C}}}{\left(\tau_{\text{C}}\right)^{2}} \tag{2}\]  
where, \(E\) is the transversal elastic modulus and \(\gamma\) is a nondimensional parameter that depends on the elasto-plastic response of the material. Different models are available in literature to estimate the value of the parameter \(\gamma\), in this work a value of \(\gamma=0.884\) has been adopted [30]. The conventional assumption is that the fracture toughness of a composite material remains constant and independent of the crack length, and mechanisms like fiber bridging, which can modify fracture propagation characteristics, are frequently neglected. However, this assumption is only valid for brittle crack propagation where the size of the process zone (\(l_{pz}\)) is negligible compared to the others characteristic dimensions, like the crack length. In practice, as the crack develops, due to the already mentioned toughening mechanisms, the formation of the process zone behind the crack tip is gradual, resulting in an increase of the fracture toughness until its value stabilize (Figure 2). Under the assumption that the cohesive stresses are independent from the profile of the crack opening and are a function of only the distance from the crack tip, an approximate expression for the resistance curve is developed in [5] using a linear idealization:  
\[G_{R}= \begin{cases}G_{T}+G_{b}\frac{\Delta d}{l_{pz}}\left(2-\frac{\Delta a}{l _{pz}}\right)&for\quad\Delta a<l_{pz}\\ G_{T}+G_{b}&for\quad\Delta a<l_{pz}\end{cases}