The fracture mechanics-based virtual crack closure technique (VCCT) is commonly used to simulate composite delamination failure. However, this technique is limited by the need for a pre-defined crack, solution convergence issues, and computational inefficiency. To overcome these limitations, an alternative delamination simulation method using a cohesive zone model (CZM) is investigated. The CZM is realized by incorporating special cohesive elements between continuum finite elements of adjacent lamina. An energy based CZM constitutive law reduces mesh dependence, but requires specification of additional failure parameters in order to adequately predict the damage process. This paper describes a modeling method for the implementation of CZM to predict delamination damage in fiber-reinforced polymer composite laminates. Procedures for calibrating cohesive element parameters with VCCT results along with computational techniques used to suppress numerical convergence issues are discussed. Laminated composite double cantilever beam (DCB) and end notch flexure (ENF) examples are analyzed using the commercial finite element code ABAQUS. Load versus deflection curves predicted using theory, VCCT, and CZM analyses are compared.