{'fraction_non_alphanumeric': 0.0217592454407068, 'fraction_numerical': 0.008894725845446677, 'mean_word_length': 5.272982587530425, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From fracture to fragmentation: discrete element modeling -- Complexity of crackling noise and fragmentation phenomena revealed by discrete element simulations

Discrete element modelling (DEM) is one of the most efficient computational approaches to the fracture processes of heterogeneous materials on mesoscopic scales. From the dynamics of single crack propagation through the statistics of crack ensembles to the rapid fragmentation of materials DEM had a substantial contribution to our understanding over the past decades. Recently, the combination of DEM with other simulation techniques like Finite Element Modelling further extended the field of applicability. In this paper we briefly review the motivations and basic idea behind the DEM approach to cohesive particulate matter and then we give an overview of on-going developments and applications of the method focusing on two fields where recent success has been achieved. We discuss current challenges of this rapidly evolving field and outline possible future perspectives and debates.


Introduction
Fracture and fragmentation of materials have been the subject to human interest for as long as we can think, mainly due to practical reasons. For centuries research on fracture was mainly driven by catastrophic failure events of engineering constructions which occurred due to the poor understanding of relevant processes [1,2,3]. Names like da Vinci, Galilei, Griffith, Weibull, Wöhler, and Inglis among others are all related to engineering solutions to fracture [2]. The nature of fracture phenomena, however, impeded systematic theoretical studies. No more than three decades ago mainstream physics slowly started to study fracture and fragmentation problems, driven by the discoveries of a young generation of researchers that made computers accessible for their research [1,3]. Based on the pioneering works of Hans J. Herrmann and others, lattice models, fuse models and meshless particle models emerged for fracture studies that -driven by the breath-taking advances of computational and algorithmic capabilities -proved to be very successful for studying fracture and fragmentation phenomena [1,2,3]. Around that time, Cundall proposed a particle method with rigid body dynamics to model fracture of frictional cohesive materials, characteristic of geotechnical applications [4]. Under the name Discrete Element Method (DEM) a group of approaches emerged that generate the motion of an assembly of particles starting from the dynamics of its constituents. The similarities of the DEM to popular methods in other fields of research like molecular dynamics [5] or smooth particle dynamics [6], lead to cross-fertilization in algorithmic development. Today DEM is a powerful tool to simulate the breaking of heterogeneous materials beyond the point of single crack growth. Various particle geometries, material response, ways to treat cohesive, repulsive behavior and of course loss of cohesion lead to a flexible tool-set of approaches. Strategies for higher order agglomeration, coupling to continuum domains or particle based fluid solvers like lattice Boltzmann extended the reach of DEM significantly [7,8,9,10,11]. Today, applications of DEM made a substantial contribution to the understanding of the mechanical response and breaking phenomena of heterogeneous materials under various types of loading conditions. Ranging from the slowly changing sub-critical loads to the highly energetic fragmentation, DEM proved to be an indispensable tool for investigations.
In this article we briefly review the motivations and basic ideas behind the DEM approach, as well as, its current extension by coupling to a continuum domain. Among the widespread applications of DEM for the fracture of heterogeneous materials we highlight two fields where recently DEM have played a decisive role to achieve major success. Finally we discuss remaining challenges of this rapidly evolving field and outline possible future perspectives.

Discrete Element Models for cohesive particulate materials
Fracture and fragmentation are difficult problems to handle numerically due to the continuous generation and evolution of crack surfaces. Classical numerical methods such as Finite Element, Finite Differences, or Boundary Elements solve partial differential equations of continuum mechanics, so that they are able to consider only a small number of discontinuities and cannot encompass the entire fracturing process. Discrete Element Modelling (DEM) is a computational approach to the deformation and failure of cohesive frictional materials which embeds materials' complexity by representing it with a set of discrete elements. The method is physically based in the sense that the elements of discretization are physical entities having e.g. mass and velocity, hence, they are called particles. The interaction of particles is defined such that the model accounts for the proper macroscopic response of the medium including both, constitutive behavior and failure mechanisms. The approach was introduced by Cundall and Strack [4] in 1979 which then initiated a rapid development of the technique and a wide variety of applications in diverse fields of engineering, physics, and geosciences [4,12,13,14,15,16].

Model construction
DEM is best suited for materials which are inherently disordered on the mesoscopic scale, i.e. they are composed of grains of various shapes with complicated cohesive coupling in between. To begin with, the model has to give a high quality representation of the microstructure of the specific material considered. To keep the problem numerically tractable, particle shapes are usually idealized by spheres in three dimensions (3D) so that the problem of discretization is reduced to the generation of a random homogeneous packing of spherical particles with a prescribed size distribution and a desired density. Two classes of generation methods can be distinguished, i.e. the dynamic and constructive methods [21,22,23,24,25,26]: Dynamic methods typically start from a random configuration of point-like particles which are then gradually blown up to the desired size reaching a dense arrangement in the domain of interest. As an alternative, particles with the required extension can be placed in a volume significantly larger than the domain of interest and then either the volume can be slowly reduced until an appropriate packing is reached or the particle system can be compactified under the action of a force field [22,13,17,19,20].
Already the packing generation involves demanding simulations of the motion of particles making dynamical methods rather time consuming. Constructive algorithms take a different strategy, namely, they are purely based on geometrical procedures to discretize the spatial domain in terms of spheres [23,26]. b) a) c) Fig. 1. Demonstration of the construction of DEMs: (a) Random homogeneous packing of spherical particles is generated by particle deposition in a rectangular container [19]. (b) Delaunay tetrahedral mesh is constructed with the particle centers and beams are introduced along the edges of tetrahedra. A close-up on the beam lattice is presented in [19,17]. (c) An early stage of the impact fragmentation of a rectangular sample generated by a projectile which hits the middle of the front side of the body. Beams are colored according to the axial strain where stretched and compressed beams have red and green colors, respectively. Colors are randomly assigned to fragments.
Efficient algorithms have been developed to fill containers of various shapes which all share the feature that they lead to packings with a low coordination number and a high porosity. The density can be increased by using the sedimentation technique [22,21,23] or by gradual refinement of the packing using tetrahedral meshes [23].
Under external load, the particle ensemble deforms and cracks emerge at highly stressed locations, the physics of which has to be captured by the dynamics of inter-particle contacts. Since the numerical representation of the deformation of particles is computationally not feasible, contact models rely on the overlap of the spherical particles and express the normal component of the contact force in terms of the overlap distance. In this so-called soft particle contact model, tangential forces and torques depend on the relative displacement of the particles since contact has been established. Realistic contact models capture dissipation, rolling and torsion resistance, as well as elasticplastic contact deformation [27,28,17,24,18].
Cohesion of the material arising due to bonding of its grains can be captured by coupling neighboring particles via elastic spring or beam elements. In the most realistic case beams in 3D account for the stretching, compression, shear, bending, and torsion of cohesive contacts [13,18,17]. Beam elements may act as bonds coupling either the surfaces [13] or the centers of mass of particles [14,15,17,19]. The geometrical features of beams, i.e. length, cross-sectional area, and moments are determined by the particle packing which leads to disorder in the bond network. The macroscopic response of the model is mainly determined by the constitutive laws and breaking criteria of beams which have to be chosen to account for the observed materials' behavior.
The primary fracture mechanism of cohesive frictional materials is the tensile and shear failure of bonds along the grain boundaries. This is captured by DEM approaches such that failure criteria of beams are formulated in terms of axial stresses and bending and torsion moments (strains) [4,12,13,14,15]. Cracks form due to the gradual removal of cohesive elements as they fulfill the failure condition during the time evolution of the system. The structural disorder of the particle packing and bond network can be complemented by strength disorder treating the parameters of the failure criterion as stochastic variables [17]. Contact forces between particles are set on when cohesion is lost to prevent the penetration of crack faces into each other. The time evolution of the particle system is followed by molecular dynamics simulations, i.e. the equation of motion of all particles is solved numerically for the translational and rotational degrees of freedom with properly set initial, boundary, and loading conditions [5,18]. The model construction is illustrated in Fig.  1 where a sedimentation algorithm was used to generate the initial particle packing ( Fig. 1(a)) [19,20]. Delaunay partitioning was carried out with the particle centers and beams were introduced between particles along the edges of tetrahedra ( Fig. 1(b)). Finally, the model was applied to investigate the impact induced breakup of a rectangular specimen ( Fig. 1(c)).

Concurrent Discrete/Finite Element coupling
Fracture in heterogeneous materials can also be understood by the flow of elastic energy from a volume into the formation of new internal surfaces [7]. Crack growth is thus a localization phenomenon with a spatially limited process zone. Physical access to the dynamic processes inside this zone can be obtained on a mesoscopic level by DEM simulation with a sufficient number of elements. However the process zone is embedded in an elastic foundation and usually the majority of particles is needed for representing the elastic domain, a job that can much more efficiently be dealt with by continuum methods. The last decade has seen an avalanche of works on different multi-scale methods for all kinds of applications and methods that can be classified to be either hierarchical or concurrent. The latter ones embrace all approaches where a fine-scale model is embedded and intimately coupled to a coarse-scale model like the example shown in Fig. 2. For a comprehensive review on the methods, we refer to [8]. As we calculate dynamic interactions in the DEM, the challenge is to couple the DEM domain to the continuum domain in a way that the interface is without spurious reflections, or it other words "mechanically transparent". Since both methods discretize time and space, they can only resolve oscillations up to a cutoff frequency w max defined by the ratio of the wave speed with respect to the minimum node distance or characteristic particle size, respectively [8]. In general the cutoff frequency of the continuum domain w c max >> w DEM max in order to benefit from a continuum approach.  Unfortunately, this results in phonon reflections at the model interface for frequencies below w c max that need to be mitigated in one way or the other.
In principle one can impose a direct edge-to-edge or master-slave coupling and damp the reflected phonons close to the model interface in the DEM domain to obtain "silent boundaries" (see Fig. 2). This master-slave coupling is a standard technique in FEM [9] and compatibility is enforced for all coupled degrees of freedom by constraining and mapping the slave DEM nodes onto the respective FEM master surface positions by the shape functions of the used elements. The forces and moments from the DEM nodes are in return added to the continuum model by standard contact procedures. Alternatively one can impose a smooth transition between models with an overlapping or bridging domain. The bridging domain method, proposed by Belytschko and Xiao [10] avoids sharp interphases by enforcing compatibility inside an overlapping domain by Lagrange multiplies. Both methods are schematized in Fig.  3. The linear scaling of relative importance of energy contributions of the different domains in the overlapping one by the blending function α assures a smooth transition between the domains.
The shape functions of the Lagrange multiplier mesh, as well as element types and geometries can differ from those of the continuum mesh. In the simplest case linear functions on triangular elements that mesh the overlapping domain are chosen with linear blending functions, but Dirac delta functions and higher order blending functions are reported to work best [8]. The reason is that these strict Lagrange multipliers enforce exact compatibility with the finite element approximation and therefore fulfill the patch test, while other types of interpolations via shape functions dont and ghost forces exist. To assure a smooth transition, also the extension of the Lagrange multiplier field should be chosen such, that several FE nodes are captured (see Fig. 4).

Rupture cascades in the discrete element model
Since the end of the '70s DEM gained widespread applications and had a substantial impact on our understanding of fracture processes of heterogeneous materials. In the following we highlight two aspects of breakdown phenomena where the heterogeneity of materials plays a crucial role and recently DEM simulations combined with the approach of statistical physics led to new understanding. First, we present how the statistics and dynamics of crackling noise emerging in slowly compressed brittle materials can be investigated in the framework of DEM, then we focus on the DEM modelling of fragmentation processes induced by energetic loading.

Crackling noise during compressive failure
Macroscopic failure of heterogeneous materials under slow external driving, i.e. under slowly increasing deformation or force, occurs as the culmination of damage accumulation [1]: at the beginning of the loading process micro-cracks nucleate at the weakest points in an uncorrelated way. Later on as the local stress fields of such defects interact, spatial correlation develops, which leads to growth of existing cracks and to an enhanced nucleation in their vicinity. The final stage of the process is dominated by the merging of cracks leading to the emergence of a macroscopic crack which spans the entire sample. However, this damage accumulation is not a "smooth" process, it proceeds in bursts of cracking events on the micro and meso scales. Such intermittent breaking avalanches generate elastic waves which can be recorded in the form of acoustic noise [1,2,3]. The acoustic emission technique is one of the most important diagnostic tools providing very valuable information about the microscopic dynamics of fracture [29]. The recent progress achieved in experimental techniques addressed the question whether crackling noise measurements could be used to forecast the imminent catastrophic failure event. The problem has a high importance for the safety assessment of engineering constructions and for the forecasting of natural catastrophes such as landslides and earthquakes [30,31,32].
Statistics of breaking bursts is usually investigated in the framework of stochastic lattice models, which are based on regular lattices of springs, beams, fibers, or fuses [1,3]. They have the advantage of allowing for a straightforward identification of breaking avalanches, however, they impose simplifications on the micro-structure of materials and on the dynamics of local breakings. Stochastic lattice models have qualitatively reproduced the integrated power law statistics of crackling noise and revealed interesting effects of the amount of disorder, friction, and range of load redistribution on the value of the exponent [3,33,54].
Both, under field conditions and in engineering applications, materials are often subject to compressive loading. Hence, the computational modelling of crackling noise under compression have a high practical importance, however, in this case lattice models of fracture face difficulties to fully capture the relevant microscopic mechanisms. To overcome this problem, recently, a DEM approach has been proposed to investigate the dynamics and statistics of rupture cascades [19,20].

Cascades of beam breaking
In DEM crackling bursts are identified as cascades of micro-fractures, i.e. correlated trails of breaking particle contacts which makes it possible to decompose the process of damage accumulation into a time series of elementary events of fracturing [19,20]. In these studies strain controlled uniaxial compression of cylindrical samples was simulated (see Fig. 5) measuring the macroscopic response of the system and the microscopic evolution of damage. A representative example of the constitutive curve σ(ε) of the system is shown in Fig. 5 is that after a beam breaks, the released stress must get redistributed in the surrounding volume. It enhances the load on neighboring beams which may induce further breakings and in turn it can even trigger an entire avalanche of breakings. In a DEM framework the breaking sequence of beams can be traced by recording the time t b i and position r b i of single breakings. In order to quantify the temporal clustering of breaking events it is assumed that consecutive breaking events are correlated if they follow each other within a correlation time t c , i.e. if the condition t b i+1 − t b i < t c is fulfilled. The value of the correlation time t c can be physically motivated, namely, it is the time needed for the stress release wave to cross the specimen.
Based on the concept of correlated breakings, bursts of local failure events can be identified. The burst size ∆ is the number of beams breaking in the avalanche which is related to the new crack surface created by the burst. It can be observed in Fig. 5(b) that during the loading process the size of bursts ∆ has strong fluctuations due to the quenched structural disorder of the material but its average has an increasing tendency towards failure. This generic behavior is in a nice qualitative agreement with the outcomes of acoustic emission measurement on heterogeneous materials [3,29,30,31,32]. In the framework of DEM, further useful quantities can be defined to characterize single crackling avalanches and the evolution of their time series: Besides the burst size ∆, the time of occurrence t and duration T are of particular interest together with the amount of energy E released by bursts. The temporal sequence of avalanches can be characterized by the waiting time t w between consecutive events.

Statistics of crackling events
The integrated statistics of the characteristic quantities, i.e. the probability distributions of the burst size ∆, energy E, and duration T , furthermore, of the waiting time t w -considering all events up to failure -proved to have a power law functional form with a stretched exponential cutoff Here x is a generic notation for ∆, E, T , and t w , while  [19,20] have an excellent agreement with the experimental findings on the statistics of acoustic bursts accompanying the compressive failure of sedimentary rocks such as sand stones [30,31,32].
Based on the detailed information DEM provides about the evolution of the crackling time series, it is also possible to investigate how the statistics of crackling events changes as the system approaches macroscopic failure.  [19,20]. This behavior is in an excellent qualitative agreement with the so-called "b-value" anomaly observed for earthquakes and in laboratory experiments on compressive fracture of rocks, i.e. the exponent b of the magnitude distribution of crackling events decreases when approaching the critical point of global failure [30,31,32].
Recent simulations also demonstrated the potential of DEM to investigate the spatial structure of damage, the gradual emergence of spatial correlation of consecutive events, the formation of the damage band due to the dominance of shear in the failure process, and even the gradual fragmentation of pieces in the damage band [19,20].

Fragmentation phenomena
Energetic loading leads to fragmentation with a multitude of cracks forming simultaneously. This leads to a rapid disintegration of solids into a large number of pieces. On a longer time scale repeated loading or shearing under a high pressure give rise to a similar outcome with fragment sizes spanning a broad range with a scale free probability distribution [34,35,36]. In Nature fragmentation of solid bodies occur on a broad range of length and time scales from the collision induced breakup of asteroids down to the degradation processes in a fault gauge. Detailed knowledge on fragmentation is required in the industry where it is exploited by technologies of mining and ore processing.
In particular such applied but also fundamental questions on fragmentation processes are most suitably answered by DEM simulations.

Universality in fragmentation
The most remarkable feature of fragmentation phenomena is that the value of the power law exponent of the size distribution of pieces shows an astonishing robustness being independent of the way of loading, of material properties, and relevant length scales. During the past decades the understanding of the observed universality has been the main driving force of fragmentation research. Experimental and numerical investigations have revealed that the universality classes of fragmentation phenomena are mainly determined by the dimensionality of the system [34,35,36,41] and by the brittle/ductile character of the mechanical response of the material [48]. For brittle materials the underlying breakup mechanisms originate from crack tip instabilities that lead to repeated crack branching-merging [41]. Combining the branching-merging scenario with the Poissonian nature of the initial nucleation of major cracks a complex functional form was proposed which describes the complete mass/size distribution of fragments p(m) including the cutoff regime, as well [41] p(m) ∼ (1 − β)m −τ exp (−m/m 0 ) + β exp (−m/m 1 ).
Here, τ denotes the exponent of the power law regime, m 0 and m 1 are characteristic fragment masses, and β controlls the contribution of the two terms of the right hand side [41]. The universality of fragment mass/size distributions is demonstrated in Fig. 7 for closed shells in 3D where shells made of three different materials were fragmented by explosion and impact against a hard wall. In the regime of small fragment masses best fit was obtained with Eq. (2) using a unique exponent τ = 1.35 ± 0.02 which defines the universality class of brittle shells. High speed imaging of shell fragmentation provided direct proof of the predicted breaking scenario [44], and additionally it revealed that not only fragment sizes but even the shape of fragments obeys scaling laws [40].

DEM simulations of fragmentation processes
Due to the difficulties of the experimental investigation of fragmentation phenomena, DEM simulations had a major contribution to the development of the field. Although, recent applications of high-speed cameras and 3D imaging have allowed for a deeper insight into the process of the rapid breakup of solids [37,38,39,40], DEM simulations still complement the experiments with very valuable information. Energetic loading like explosion or impact generate a large number of simultaneously growing cracks which interact with each other in a complicated way. DEM has the capabilities to handle this high degree of complexity making possible a realistic treatment of fragmentation processes.
Detailed studies with DEM in various embedding dimensions revealed a transition from damage to fragmentation [45] at a critical imparted energy already for two-dimensional systems. The existence of the damage-fragmentation critical point has been confirmed by further DEM simulations [46,47]  In a similar impact simulation of a plastically deforming sphere at low impact velocities, a single crack occurs in the middle and a large permanent deformation remains around the impact site [48]. (c) Impact experiments of plastic balls revealed a similar breakup mechanism in a good quantitative agreement with DEM simulations [48].
types of fragmentation processes, and it was also reproduced by experiments [42,43]. The result implies that universality of fragment size distributions is due to the underlying continuous phase transition. The entire richness of fragmentation mechanisms however could only be resolved by full 3D systems once a significant particle number could be considered [17]. Contrary to the simple branching-merging scenario it became clear that there exist different mechanisms which get activated as the imparted energy increases and their interaction determines the final breaking scenario [17]. Fragmentation processes of plastically deforming materials show an even higher complexity: power law size distribution of fragments has been confirmed by experiments, however, with an exponent significantly lower than for brittle materials. DEM simulations clarified that shear induced breaking is responsible for the emergence of the novel universality class which makes plastic fragmentation similar to the one of liquid droplets [48]. The effect of material microstructures on the outcome of the fragmentation was studied by mapping material micro-structures of multi-phase materials and composites onto the DEM systems [49]. Surprisingly the size distribution exponent is rather robust with respect to such issues, strengthening the universal character of fragmentation. Representative examples of DEM simulations of fragmentation processes are presented in Fig. 8 both for brittle and ductile materials. The figure also demonstrates the agreement of simulations with experiments.

Conclusions and future challenges
In the present paper we briefly reviewed the basic ideas behind DEM for heterogeneous materials and highlighted two fields where this modelling approach played a decisive role to reach recent success. Profiting from the increase of computer power and the success of hybridization of modelling approaches simulation studies of fracture and fragmentation phenomena can help to resolve current debates of the field and to reach new challenges. When studying statistical features of the fracture of heterogeneous materials such as the size effect of macroscopic strength and crackling noise generated by avalanches of microfractures, stochastic lattice models have been successfully applied under tensile loading conditions [54]. However, under compressive loading they usually have difficulties to account for all relevant mechanisms. We have demonstrated that DEM offers an adequate modelling framework for crackling phenomena under compression reproducing all observed scaling laws of rupture cascades obtained by field measurements, as well as, in laboratory experiments [19,20]. The results imply that DEM has a high potential to understand the emergence of catastrophic failure in porous granular media challenging Earth sciences and engineering. In natural catastrophes such as landslides and earthquakes, the available data are often incomplete and provide only a limited insight into the complexity of processes that lead to failure. The main contribution of DEM is that it can capture all relevant processes down to the length scale of single grains, and hence, it can reveal mechanisms hidden for experimental approaches. With these capabilities for the investigation of the statistics and dynamics of rupture cascades DEM may give rise to a breakthrough in developing predictive models of catastrophic failures and even earthquakes in the near future.
Research on fragmentation faces similar challenges. Recent experiments on impact induced fragmentation of one-and two-dimensional objects revealed that the power law exponent τ of the fragment size/mass distribution increases with the imparted energy [50]. DEM simulations performed in two dimensions confirmed this finding and yielded a logarithmic dependence of τ on the energy [51,52]. The results are remarkable because on the one hand they question universality and the underlying phase transition picture of the damage-fragmentation transition, and on the other hand, they have relevance for industrial applications in mining and ore processing, as well. However, recent DEM studies on the breakup of spherical bodies due to impact against a hard wall demonstrated that the apparent increase of the exponent can be removed by rescaling the mass/size distributions with the average fragment mass [46]. Finite size scaling proved to be indispensable to correctly determine characteristic exponents of fragmentation phenomena which again shows the importance of large system sizes and calls for further investigations to settle the problem. Both experimental and theoretical investigations have shown that the dimensionality of the breakup process, especially the interplay of the dimensionality of the object and of the embedding space, plays a crucial role in the selection of the dominating mechanism of dynamic cracking and fragment formation [44,41,53]. This addresses the opportunity that in certain cases universality can be violated and the energy dependence can be understood through the gradual activation of different breakup mechanism and the mixing of them as the imparted energy is varied.
Comminution of solid particles in ore processing is one of the most important applications of fragmentation inducing a huge cost for the industry. Improving the efficiency of size reduction in comminution machinery and the separation of fragments according to their size are still important driving forces of fragmentation research greatly profiting from DEM.
Advancement of measuring technologies has made it possibly to go beyond the analysis of the mass/size distribution of fragments in the final state of the breakup process. There is an increasing amount of information available on the velocity of pieces, as well. It is a great challenge for theoretical investigations to understand what determines the functional form of the velocity distribution of fragments, and whether the mass and velocity of fragments are correlated. Beyond their scientific importance both problems have also practical relevance: on orbit fragmentation events are the main source of space debris where the velocity of debris pieces and the presence of mass-velocity correlation are crucial for estimating the risk of damaging collisions with satellites.
Today particle models for fracture and fragmentation are at the verge of becoming significant tools for simulating industrial processes. The dilemma of either using a large number of spherical particles or a significantly smaller number of aggregated or polygonal particles for discretization is slowly diluted by the development of computer hardware and should vanish within the next decade. Additionally, the incorporation of DEM into FEM workbenches will bring these methods to a wider community of applied users. With every new release of FEM simulation suites, software companies extend functionality, recently to incorporate particle methods. Even though still the simplest methods are implemented, soon we might see advanced DEM embedded in FEM code with robust concurrent coupling. As the continuum and discrete worlds continuously merge inside commercial software packages users are increasingly liberated from technicalities of discretization and implementation issues. Till today a major problem of engineering design is the fracture size effect, i.e. the fracture strength of heterogeneous materials decreases with increasing sample size [54]. The size effect is the consequence of the complex interplay of the long range elastic interaction of material elements and of the inherent disorder. Since DEM naturally accounts for both, with the increase of computer power and algorithmic developments DEM is expected to have important contribution to tackle fracture size effect.