A constitutive model for an interface finite element is proposed to enable simulation of delamination in composite materials with R-curve effects. The constitutive model is formulated in the framework of cohesive zone modelling (CZM). In essence, a multilinear CZ law with an arbitrary number of line segments is developed. The CZ law seeks to enable constitutive modelling of failure mechanisms on multiple scales within the fracture process zone and reduce conventional a priori assumptions regarding the shape of the CZ law. The CZ law relies on damage mechanics, an equivalent one-dimensional formulation, and criteria for mode interactions to simulate delamination under mixed-mode loading. Special emphasis is put on the derivation of interpolation formulas and a constitutive tangent stiffness tensor for the multilinear formulation. The constitutive model is implemented in the commercial FE program ANSYS Mechanical, for implicit finite element analysis (FEA), using user-programmable features. The implementation is verified through single interface element numerical studies, and its applicability is demonstrated by simulating an experiment of quasi-static delamination showing large-scale fiber bridging in pure mode I DCB glass-fiber epoxy specimens. Experimental measurements and simulation outputs using the novel cohesive element is compared to those of the conventional bi- and trilinear CZ laws.