Debonding of the adjacent layers in composite materials, or delamination, is one of the most common and important damages which could lead to an abrupt decrease in the strength and stiffness of the structure. Delamination is driven by either service life loads or manufacturing processes and has been widely investigated by the Cohesive Zone Models (CZMs) and interfacial elements considering all bulk and interface layers engaged [1]. Since such an approach could be inefficient specially at large number of layers, a method based on the homogenization theory through the Rule of Mixtures (RoM) is presented in this study, where delamination damage is introduced at the composite level in a way that the bulk response of the composite is properly simulated. To this end, the parallel mixing theory [2] has been employed to reproduce the bulk response of the composite based on an iso-strain assumption. Potential delamination areas are then identified once interfacial shear and normal stresses reached the threshold values, and the delamination damage develops based on a stress formulation [3]. To assess the implementation performance, several simulations were conducted in different mode mix ratios through Double Cantilever Beam (DCB), End notch flexure (ENF) and Mixed mode bending (MMB) tests. Arc length strategy was employed to run the cases according to the snap back phenomenon whilst element deactivation was utilized to remove the element after complete failure. Results showed that the adopted method could be considered as a promising way to incorporate delamination damage effects within the conventional finite element methods.