A cohesive zone model implemented in an augmented Lagrangian functional is used for simulation of meso-scale fracture problems in this paper. The method originally developed by Lorentz is first presented in a rigorous variational framework. The equivalence between the stationary point of the one field problem and the saddle point of mixed formulation is proved by solving the double inequality of the mixed functional. An adaptation to simulate fracture phenomenon in the meso-scale via mesh modification is also presented as an algorithm to insert zero-thickness interface elements based on Lagrange multipliers, boarding the non-trivial task of the field interpolation for different crack paths (plain and tortuous). A suitable tool to study the matrix fracture and debonding phenomena in composites with strongly different component stiffnesses that avoids ill-conditioning matrices associated with intrinsic cohesive zone models is obtained. The method stability is discussed using a simple patch test. Some numerical applications to fracture problems taking into account the meso-structure and, particularly, the study of transverse failure of longitudinal fiber reinforced epoxy and the fracture in concrete specimens are included in the paper. Comparing the numerical results with the experimental results obtained by other researchers, the paper introduces a discussion about the influence of coarse-aggregate volume in meso-scale fracture mechanism in concrete L-shape specimens.