Design and maintenance of concrete structures require efficient methods for predicting quasi-brittle fracture while incorporating mesoscale-induced structural uncertainties. Toward this aim, a phase field cohesive zone model combined with the cell-based smoothed finite element method (CSFEM-PFCZM) is developed in the study, with mesoscale heterogeneity of concrete characterised by Weibull random fields. An automatic quadtree decomposition technique is presented to adaptively refine meshes around mesoscale areas, in which the hanging-node problem is tackled by treating hierarchical quadtree elements as CSFEM polygons. The crack driving forces are evaluated as history variables at the integration points of CSFEM subcells, and only the boundary geometry and Wachspress shape functions are needed to calculate strain matrices and stiffness matrices without the Jacobian inversion problems. The developed method is validated by extensive Monte Carlo simulations of typical nonlinear fracture benchmarks of concrete structures, demonstrating that the well-captured variability of crack paths and load–displacement curves can supplement the limited statistical results obtained from a small number of experiments. It is also found that larger correlation length and higher variance underline higher heterogeneity that leads to more dispersed responses. The developed method holds promise for the efficient evaluation of structural reliability taking into account stochastic mesoscale effects, due to the fast generation of random field samples in large quantities and flexible simulation of complicated fracture through the CSFEM-PFCZM.