The mechanical performance of concrete is strongly influenced by the geometry and properties of its components (namely aggregate, mortar, and Interfacial Transitional Zone (ITZ)) from the mesoscale viewpoint, and analyzing the material at that level should be a powerful tool for understanding macroscopic behavior. In this paper, a simple and highly efficient method is proposed for constructing realistic mesostructures of concrete. A shrinking process based on 3D Voronoi tessellation was employed to generate aggregates with random polyhedron and grading size, and reversely, an extending procedure was applied for ITZ generation. 3D mesoscale numerical simulation was conducted under a quasi-static load using an implicit solver which demonstrated the good robustness and feasibility of the presented model. The simulated results resembled favorably the corresponding experiments both in stress–strain curves and failure modes. Damage evolution analysis showed that the ITZ phase has profound influence on the damage behavior of concrete as damage initially develops from here and propagates to mortar. In addition, it was found that tensile damage is the principal factor of mortar failure while compressive damage is the principal factor of ITZ failure under compression.