## 1 Introduction
The numerical approaches adopted to simulate delamination propagation in Finite Element (FE) analysis typically use the value of the fracture toughness measured at the onset of the delamination, underestimating the load-carrying capability of the structure and resulting in conservative designs, as the increasing resistance of the delamination, the R-curve effect, is ignored [16]. Delamination propagation can be numerically simulated in FE analysis using approaches based on fracture mechanics, such as the Virtual Crack Closure Technique (VCCT) [17; 18], or damage mechanics, such as the Cohesive Zone Modelling (CZM) [19; 20; 21]. In the last decades, both these approaches have been modified to take into account the R-curve effects. Several researchers have demonstrated the possibility of using VCCT-based methodologies by modifying the fracture toughness value used in the propagation criteria following the experimentally determined curves of the fracture toughness as function of the crack opening [22; 23; 24]. These approaches can be easily incorporated in a FE analysis, however their applicability is limited to scenarios where the location of the initial delamination front and the direction of propagation are known a priori. On the other hand, conventional cohesive zone model (CZM) using a bilinear traction-separation law has been developed and extensively utilized in numerical simulations to study delamination initiation and growth, without considering fiber bridging. Subsequently, various modifications to the traction-separation law have been proposed to consider the additional fiber bridging effect based on the R curve obtained experimentally. [PERSON] et al. [25; 26; 27] proposed a tri-linear law, obtained from the combination of two bilinear cohesive laws to simulate the toughening mechanisms during the propagation of the delamination. This approach has been adopted with different modifications by other researchers [28; 29] demonstrating its effectiveness in considering the effect of the R-curve in standard specimens. Although the use of cohesive elements is more computationally expensive compared to VCCT, they enable the simulation of delamination initiation without the need to define an initial crack tip. Approaches based on the cohesive zone model to represent the R-curve effect have been mainly applied to standard specimens, such as the Double Cantilever Beam (DCB). The objective of this research is to demonstrate the applicability of this numerical technique to a more complex specimen. In this work, the effect of R-curve on the delamination propagation in skin-doubler specimens, representative of the stringer flange of stiffened composite panel, subjected to three-point bending test, is numerical investigated using a superposition of two bilinear cohesive law. The specimens were manufactured with the same material exhibiting two different fracture behavior: one with limited R-curve effect, indicated as \"Low Bridging\", and the other with a significant increase in the fracture toughness as the delamination propagate, indicated as \"High Bridging\". The numerical results are compared with the experimental data in terms of ultrasonic C-scan images and load vs. applied displacement curves obtained in [23]. The comparisons validate the effectiveness of the adopted approach in replicating the increase in strength due to the R-curve effect observed experimentally and in predicting the evolution of the delaminated area as a function of the applied displacement. In Section 2, the theoretical background of the cohesive zone model employed to represent the R-curve effect is introduced. In Section 3, the geometrical characteristics of the skin-doubler specimens are described together with the details of the FE model, while in Section 4, the calibration of the cohesive zone model parameters is performed for the materials under consideration. Finally, in Section 5, the numerical results are presented and compared with the experimental data.