By utilizing numerical models and simulation, insights about the fracture process of brittle heterogeneous materials can be gained without the need for expensive, difficult, or even impossible, experiments. Brittle and heterogeneous materials like rocks usually exhibit a large spread of experimental data and there is a need for a stochastic model that can mimic this behaviour. In this work, a new numerical approach, based on the Bonded Discrete Element Method, for modelling of heterogeneous brittle materials is proposed and evaluated. The material properties are introduced into the model via two main inputs. Firstly, the grains are constructed as ellipsoidal subsets of spherical discrete elements. The sizes and shapes of these ellipsoidal subsets are randomized, which introduces a grain shape heterogeneity Secondly, the micromechanical parameters of the constituent particles of the grains are given by the Weibull distribution. The model was applied to the Brazilian Disc Test, where the crack initiation, propagation, coalescence and branching could be investigated for different sets of grain cement strengths and sample heterogeneities. The crack initiation and propagation was found to be highly dependent on the level of heterogeneity and cement strength. Specifically, the amount of cracks initiating from the loading contact was found to be more prevalent for cases with higher cement strength and lower heterogeneity, while a more severe zigzag shaped crack pattern was found for the cases with lower cement strength and higher heterogeneity. Generally, the proposed model was found to be able to capture typical phenomena associated with brittle heterogeneous materials, e.g. the unpredictability of the strength in tension and crack properties.