In the present work, a numerical model based on the cohesive zone modeling (CZM) approach has been developed to simulate mixed-mode fracture of co-consolidated low melt polyaryletherketone thermoplastic laminates by considering fiber bridging. A modified traction separation law of a tri-linear form has been developed by superimposing the bi-linear behaviors of the matrix and fibers. Initially, the data from mode I (DCB) and mode II (ENF) fracture toughness tests were used to construct the R-curves of the joints in the opening and sliding directions. The constructed curves were incorporated into the numerical models employing a user-defined material subroutine developed in the LS-Dyna finite element (FE) code. A numerical method was used to extract the fiber bridging law directly from the simulation results, thus eliminating the need for the continuous monitoring of crack opening displacement during testing. The final cohesive model was implemented via two identical FE models to simulate the fracture of a Single-Lap-Shear specimen, in which a considerable amount of fiber bridging was observed on the fracture area. The numerical results showed that the developed model presented improved accuracy in comparison to the CZM with the bi-linear traction–separation law (T–SL) in terms of the predicted strength of the joint.