Generation and packing algorithms are developed to create models of mesoscale heterogeneous concrete with randomly distributed elliptical/polygonal aggregates and circular/elliptical voids in two dimensions (2D) or ellipsoidal/polyhedral aggregates and spherical/ellipsoidal voids in three dimensions (3D). The generation process is based on the Monte Carlo simulation method wherein the aggregates and voids are generated from prescribed distributions of their size, shape, and volume fraction. A combined numerical-statistical method is proposed to investigate damage and failure of mesoscale heterogeneous concrete: the geometrical models are first generated and meshed automatically, simulated by using cohesive zone model, and then results are statistically analysed. Zero-thickness cohesive elements with different traction-separation laws within the mortar, within the aggregates, and at the interfaces between these phases are preinserted inside solid element meshes to represent potential cracks. The proposed methodology provides an effective and efficient tool for damage and failure analysis of mesoscale heterogeneous concrete, and a comprehensive study was conducted for both 2D and 3D concrete in this paper.