Unlike most cohesive zone models (CZMs), the Park–Paulino–Roesler (PPR) cohesive fracture model has the inherent ability to control the softening shape of the traction–separation relationships, which makes it suitable to model fracture for a wide variety of materials. Like other CZMs, the PPR model is well-suited for problems where the crack path is known a priori but its implementation may become complex in situations where the crack path is not known beforehand. To overcome this limitation, we recast the PPR model within the framework of the phase field method, which enables the modeling of crack propagation problems with complex crack topologies using a straightforward multi-field finite element implementation. We use constitutive functions (i.e., crack geometric function and degradation function) consistent with the PPR model for the case of mode-I fracture, such that equivalent traction–separation relationships from the phase-field model approximate those from the original PPR model. Our choice of geometric and degradation functions is based upon those proposed by Wu (2017). We present several numerical examples to demonstrate the ability of the model to capture fracture of problems with different materials, geometries, and boundary conditions. Also, we show that the results from our model converge to those obtained with the original PPR model for problems where the crack path is known a priori. Being crucial in engineering design, we finally show that the model can capture size and boundary effects with satisfactory accuracy.