Discrete element method (DEM) was used to model progressive delamination of fiber reinforced composite laminates. The anisotropic composite plies were constructed through a hexagonal packing of particle elements. Contacts between the particles were represented by parallel bonds with the verified normal and shear elastic properties. The ply interface was characterized by a contact softening model with a bilinear elastic behavior which is similar to the cohesive zone model in the continuum mechanics. DCB. ELS and FRMM tests were simulated by the DEM model to assess its capability of modeling mode I, mode II and mix mode fracture of delamination, respectively. Good agreements were observed between the DEM and existing numerical and experimental results of loading curves, which confirmed that the DEM model can be used to simulate initiation and propagation of composite delamination, with more insights into microscopic material behavior, such as damage extension and plastic zone.