Strength of industrial sandstones modelled with the Discrete Element Method


 Sandstone is one of the most popular building materials since the earliest times. It has various textures and colours as well as good technical parameters. Sandstones, having such wide applications, are subjected to various external factors during the period of use. So, it is of utmost importance to have a good knowledge of their strength parameters. We employed a numerical method called Discrete Element Method to examine in a non-invasive manner the mechanical strength of industrial sandstones, that are commonly used as broken stones in road construction, cladding material, paving stones, pavement tiles and so on. Various mechanical external factors were considered, such as breaking, compressional and abrasion forces or impact by external objects and vibrations. Fragmentation of the considered sandstones under compressional regime was a source of knowledge about energy storage inside the material and energy release, as well as appearance of fractures inside the matter and final sandstone fragmentation into crumbs.


Material failure
Strength of solid bodies is determined by the structure of the given material and the way of applying external forces leading eventually to its destruction and fragmentation (Jaeger et al. 2007; Carmona et al. 2008). Typically, fragmentation process is described at a macroscopic level (typical for engineering approaches), where fragmentation of a material can be understood as the development of a single crack in the material that finally leads to its splitting into parts (Griffith 1921;Orowan 1949). However, when it comes to rocks modelling, grainy structure of the material should also be taken into account. In such case, a microscopic approach is quite attractive because it operates at the level of grains forming the material. In the present study, a microscopic approach has been applied in an attempt to numerically model selected parameters of different sandstones. Sandstone is a type of sedimentary rock composed mainly of sand-sized (0.06 -2 mm) particles or fragments (Łukaszewski 2003, 2013; Klemm and Wiggins 2016) that relatively easily crumble and fragment due to the activity of external forces. Numerical approach was aimed at recreating samples with the same dimensions and particle sizes as real laboratory equivalents.

Numerical modelling of materials
Numerical modelling of materials generally includes two branches -one of them assumes continuity of the medium and the other its discrete structure. Both have their advantages and disadvantages (Courtney 1990). Classical engineering approach to the fracturing problem is rooted in a continuum mechanics and is mainly based on the computer technique called Finite Element Method (Munjiza 2004). However, in this approach, the discrete nature of the granular medium and the relative displacements and rotation of particles inside the material are not considered (Onate and Rojek, 2004). Particle-based numerical methods offer an alternative approach, among which the Discrete Element Method (DEM) is one of the most popular (Egholm 2007;Rojek 2007). DEM is well suitable for following the evolution of fractures in brittle materials (Rojek 2007). In DEM, the material is represented as an assembly of interacting particles -discrete elements (Cundall and Strack 1979;Cundall 1971Cundall , 1974Cundall , 1978. When two particles overlap, they start to interact and the magnitude of this interaction depends on the contact forces by the force-displacement relation (Kazerani 2013). Each discrete element can move and rotate, and contacts between elements are detected in every step of the simulation (Fraige and Langston 2004;Rojek 2007).
The algorithm of DEM ( Fig. 1) can be divided into two main parts: the first one is related to the creation of a contact model and the calculation of forces acting on the elements, and in the second part, the second Newton's law of dynamics is applied to each element to calculate the changes in position and velocity that are a result of acting forces (Cundall and Strack 1979;Cundall 1971Cundall , 1974. To maintain the numerical stability of the simulation, it is necessary to use small time steps (O'Sullivan and Bray 2004; Potyondy and Cundall 2004), so that the elements do not move too much in a time interval (element movement is small enough to affect only the immediate surroundings of the element). The DEM is computationally expensive. Even now, modelling largescale problems is extremely time consuming and requires numerous simplifications. Presented simulations were carried out with the help of the supercomputer Cray XC40 'Okeanos' of the Interdisciplinary Centre for Mathematical and Computational Modelling, Warsaw University. All the simulations were performed with open-source DEM software ESyS-Particle (Abe et al. 2014) designed especially to study the physics of rocks and the non-linear dynamics of earthquakes and to address the computational limits of existing DEM software.

Description of the model
ESyS-Particle provides different types of interactions (bonds) between particles. One of them -called BrittleBeamPrms in software notation -was designed especially for simulating brittle fracturing (Wang 2009; Wang et al. 2006). The BrittleBeamPrms interaction involves all six degrees of freedom of each interacting particle (Li 2012). Particle motion can be decomposed into two independent parts -translational motion of the centre of the mass and rotation about the centre of the mass (Fig.  2). The translational motion of the centre of the mass is governed by the Newtonian equations and integrated using a conventional Molecular Dynamic scheme (Mora et al. 1993(Mora et al. , 1994. The particle rotation is governed by the Euler's equations. Translational and rotational degrees of freedom can be included in particle-pair interactions. It means that two bonded particles interact by normal and shear forces and bending and twisting moments (Fig. 3) causing the particles to rotate relative to each other when they are in frictional contact. Particles are able to rotate  Particles a and b of radii R a and R b are at a distance D with overlap C. P a and P b are the points of intersection of the line connecting the disc centres with the boundaries of the discs; d and t are unit vectors. Derivatives of r and θ stand linear damping and rotational damping, respectively. when they are in frictional contact because rotational frictional interactions introduce a torque to both particles. Simplifying, rotational bonds are similar to cylindrical elastic beam connecting two particles. Such bond has a radius equal to the mean of the radii of the bonded particles. Equilibrium length of the bond is equal to the sum of the radii of those particles (Abe et al. 2014).
Assuming small deformations, the relationship between forces or torques and relative displacements between two bonded particles (Eq. 1) could be written in the following form: where F represents forces, M torques, K rigidity constant, ∆r and ∆s are relative translational displacement and ∆a denotes relative angular displacement. The meaning of the symbols r, s, t and b is as follows: r -radial component, s -shearing component, t -twisting component and bbending component. In the isotropic case K s = K s1 = K s2 and K b = K b1 = K b2 . Therefore only four rigidity parameters are required. The elasticity of bonds (under the above assumptions and according to the linear elastic beam theory) can be defined as a function of two parameters -the bond Young's modulus and the bond Poisson's ratio (Weatherley et al. 2010). Bond between the particles breaks according to a Mohr-Coulomb failure criterion. It happens when the shear stress inside the bond is bigger than its shear strength (τ) according to the formula τ = c + σ N tan(φ f ), where c b is the cohesion of the bond for zero normal stress (σ N ) and φ f is the internal angle of friction of the bond. Therefore, there are four bond parameters in the model to calibrate: Young's modulus E b , Poisson's ratio v b , cohesion c b and internal friction angle φ b .
Such a material for small external loading can be treated as an ideal elastic body. When loads are larger (when some bonds start to break and particles are able to displace from their initial positions), the material behaves with plasticity. Large stress finally leads to initiation and development of cracks. As a result of fragmentation, particles without bonds appear. Particles can be subjected to bonded interactions or frictional (unbonded) interactions; therefore, exclusion is introduced between these two interaction groups. Between unbonded particles, frictional interaction is applied and a broken bond is identical with a fracture surface (Abe et al. 2014).
Frictional behaviour at a microscopic scale is achieved due to a Coulomb-type friction law. The Coulomb-type friction law and the tangential force are responsible for the beginning of sliding at specific contact. This friction law also introduces static and dynamic frictional behaviours at the microscopic scale with a linear function. The magnitude of the shear force is checked against the maximum possible value (F s ) max = F n tan(φ μp ), where φ μp is the smallest of the interparticle friction angles of the two discs ( Fig. 2) in contact (Ferdowsi 2014). Friction coefficient between particles can be either static or dynamic.
Acoustic emissions appearing during fracturing are attenuated (translational and rotational oscillations) due to implemented body forces. Such so-called artificial viscosity is used to avoid unphysical accumulation of kinetic energy of the particles. However, damping should have little effect on the elastic response of the simulated material; therefore, the viscosity coefficients are small (but sufficient to damp unwanted oscillations).
BrittleBeamPrms assumes a Mohr-Coulomb failure criterion governing the maximal shear stress within a beam interaction, as a function of the normal (tensile/ compressive) stress. The formula for calculating the shear and normal stresses involves the forces and moments within the beam interaction (Eq. 2) accordingly: where τ is a shear stress, σ is a normal stress, F τ is a shear force, F σ is a normal force, A is the cross-sectional area of the beam, R is the radius of the beam, I is the bending moment of inertia and J is the twisting moment of inertia. These equations result from linear elastic beam theory and are typically employed in bonded particle DEM models involving rotational degrees of freedom (Weatherley et al. 2010). The Hoek-Brown criterion formula recreates nonlinear failure envelope. It was observed that after a series of triaxial tests on the core samples over a wide range of confining pressures, the failure envelope was found to deviate from a linear Mohr-Coulomb failure envelope for higher confining pressures. The Hoek-Brown criterion is an empirical approach for this phenomenon. When considering numerical model, the macroscopic mechanical response of the assembly of the particles is only partially determined by the choice of interaction law. Such interaction laws (like BrittleBeamPrms) are applied at a microscale, the scale of individual interactions between DEM spheres. A considerable role in the macroscopic response is also played by the network of bonded interactions. Different macroscopic responses will be obtained for different arrangements of DEM spheres or differing size distributions. The Hoek-Brown failure criterion was created to mimic the macroscopic response of rocks. It is not clear if such a criterion should be employed at a micro-scale, between two bonded DEM spheres, because there is no guarantee that it will result in a non-linear, macroscopic failure envelope. DEM investigations suggest that it will not be the case in general (Weatherley et al. 2010). A tensile cut-off can be included as an option in BrittleBeam interactions. However, it is not yet found to have any appreciable macroscopic influence on the mechanical properties of the bonded particle assembly. In the interest of keeping the model simple, it is tend not to employ such a cut-off at the scale of individual beam interactions.
Original foundations of DEM come, first of all, from Cundall (1971Cundall ( , 1974 and Cundall and Strack (1979) and further explanations can be found in Rojek (2007).

Model validating
A standard procedure in DEM modelling is model calibrating. Calibration using DEM simulations is actually an iterative process of adjusting input parameters such that the macroscopic results of simulations and experiments are equivalent. Such trial and error empirical procedures are very time consuming and unpractical and, despite some new attempts and propositions (Do et al. 2017), this is still the only most sensible procedure. In this paper laboratory geomechanical testing was used to determine the strength and deformation properties of sandstones, which were later used during the calibration process. Four types of sandstones were the object of study, with different lithology and coming from different regions of Poland (Fig. 4).
Monoliths of joint sandstone of the upper Cretaceous were collected from the active quarry in Radków in the Table Mountains (Sudetes), which were medium and coarse sandstones with silica-clay cement and yellow colour. Monoliths of an elliptical, laminated sandstone of the lower Triassic (Bunter sandstone) were collected from the active quarry in Tumlin in the Świętokrzyskie Mountains. These were medium-grained sandstones with a silica-clay-iron cement and red colour. Monoliths of Godula sandstone of Cretaceous (Cenoman) were taken from an active quarry in Brenna in the Silesian Beskid (in the Carpathians). These were coarse and mediumgrained flysch sandstones with carbonate-chlorite cement and grey colour. Monoliths of Krosno sandstone of Paleogen age (Eocen/Paleocen) were collected from the active quarry in Mucharz in the Little Beskid Mountains (in the Carpathians). These were medium-grained flysch sandstones with carbonate-silty cement and grey colour. In further investigations sandstones will be named after their place of origin: Brenna, Mucharz, Radków and Tumlin. Such materials are commonly used in industrial applications. Sandstones from Brenna are used as crushed stones in road construction and as cladding material. Sandstones from Mucharz are used as cladding material, block stone and paving stones. Sandstones from Radków are used as cladding material. Sandstones from Tumlin are used as cladding material, building elements and pavement tiles. Therefore , it is of special interest to find out how they react under different conditions. The scheme of the assay can be divided into two parts. Firstly, samples of sandstones were tested in a laboratory. Uniaxial compression strength test was employed to obtain sandstone characteristics. The obtained results were the basis for calibrating in further DEM modelling. Outcomes of laboratory and numerical uniaxial compression strength were applied to tune DEM model parameters. In the end, all calibrated numerical sandstones were subjected to different strength tests.

Laboratory investigation of sandstones under uniaxial compression strength
Uniaxial compression strength test is a common laboratory test in which a cylindrical sample is compressed from the top and bottom, and it provides a simple and effective way to characterise a material's response to loading (Hertzberg 1976;Łukaszewski 2003;Rojek 2007).
Specimens from drilled cores are prepared by cutting them to the specified dimensions, where the recommended ratio between the height and diameter of the sample is between 2 and 3 (Rojek 2007). By subjecting a sample to a controlled tensile or compressive displacement along a single axis, the change in dimensions and resulting load can be recorded to calculate the stress-strain profile. From the obtained curve (so-called stress-strain curve), elastic and plastic material properties can then be described (Rojek 2007). The following relations are used to calculate various quantities in this type of experiment. Strain can be calculated as ε e = ΔL/L 0 where ΔL is the measured displacement and L 0 is the initial sample length along a single axis. Stress is defined as σ = F/A, where F is the applied load and A is the initial cross-sectional area of the sample normal to the loading direction. Material subjected to loading initially behaves in a linear elastic manner, that is, stress and strain are linearly related. The slope of the strain-stress curve within the linear elastic regime is equal to Young's modulus, which is the ratio of the stress to axial strain E = σ/ε z (Sochor 1998). Poisson's ratio is the ratio of the relative transverse strain normal to the applied load ε xy and the relative axial strain ε z in the direction of the applied load. Poisson's ratio can be expressed as ν = ε xy /ε z .
Uniaxial compressive strength (UCS) tests were performed according to the Polish Standard PN-EN 1926 (2007) test procedure 'Natural stone test methods -Determination of uniaxial compressive strength'. Twelve samples of sandstones were tested. Drilling core sections with a diameter of approx. 38 mm were cut to a lengthto-diameter ratio equal to 2. The strength tests were performed on the samples after they were dried at 70°C to a constant weight and they reached a thermal equilibrium at 20°C for 24 hours. The strength tests were carried out in the MTS-815 rigid press produced by MTS System Corporation.
During the test, the load was applied continuously, with a constant stress rate of 1 MPa/s. For samples with a diameter of about 38 mm, the load rate was 65 kN/min. Load changes (F), axial strain (ε z ) and circumferential strain (ε xy ) were recorded. The value of volume strain (ε v ) was determined on the basis of the values of the axial strain and the circumferential strain: ε v = ε z + 2ε xy . The value of the maximum load at failure (F max ) was determined to the nearest 1 kN. Next, the value of the UCS was calculated as follows: where F max is the maximum load at failure (kN) and A is the cross-sectional area (mm 2 ).

Fig. 5: Schematic representation of the uniaxial compression test. Sample is compressed from the top and bo
On the left -simplified scheme of the test, on the right -laboratory experiment.   During the tests, experimental deformation stressstrain curves were plotted. They were used as the basis for determining the elasticity parameters according to the European Standards (Eurocode 7): the Young's modulus (E) was defined as the average modulus of the relationship between the stress σ and the axial strain ε z of the straight line portion of the stress-strain curve. Large differences were obtained for UCS. The values ranged between 55 and 153 MPa. Again, the differentiation depended on the type of sandstone. were obtained for medium-grained sandstone from Tumlin and medium-grained Krosno sandstone from Mucharz, respectively. The medium values of Young's modulus E (22.5-24.9 GPa) were obtained for coarse-grained joint sandstone from Radków. Full set of obtained parameters is presented in Table 1.
During deformation of the examined sandstones, the process of dilation was observed, that is, a relative increase in the volume of rock samples under loading. The relative increase in volume in relation to the expected elastic changes, that is, the threshold of relative dilatation (microdilation), was recorded from 40% UCS in the sandstones from Radków to 55% UCS in the sandstones from Mucharz. On the other hand, the proper increase in volume, that is, the threshold of absolute dilatation (macrodilation), from which unstable cracking already begins, was recorded from 50% UCS in the sandstones from Radków and Tumlin to 76% UCS in the sandstones from Mucharz.

Numerical investigation of sandstones under uniaxial compression strength conditions
The next step, after laboratory measurements, was DEM modelling. Four cylindrical models were created for all four types of sandstones. Every cylinder had a height of 75 mm and a radius of 18.75 mm -the same values as in case of samples used in laboratory. The presented models were constructed of 82,266 particles jointed by 297,847 bonds (Brenna), 339,324 particles jointed by 821,488 bonds (Mucharz), 150,653 particles connected by 500,024 bonds (Radków) and 1,155,851 particles jointed by 3,831,657 bonds (Tumlin). The number of particles was not an arbitrary choice but was determined by the particle packing algorithm and depended on the cylinder geometry and particle size. Applied geometric algorithm (implemented in the software) places elements according to their geometric relationships, often using pre-prepared grids. The biggest particles are firstly located in the model domain as the seeds. Then, subsequent particles are inserted in the remaining spaces. This process is completed when the model domain is filled up. The full procedure can be found in Mora and Place (2002). Therefore, the number of particles presented above is not a predetermined parameter, but a result of sample geometry. At the beginning, all particles were connected by bonds; therefore, at the initial stage, there were no microcracks inside. DEM simulations are time consuming (Klejment 2020). Computational time of a single test varied from almost 5 hours (using 93 processors of the supercomputer for a model with the smallest number of particles -Brenna), up to 48 hours (using 320 processors of the supercomputer for a model with the highest number of particles -Tumlin). Bulk density of the particles was set as 2429 kg/m 3 (Brenna), 2655 kg/m 3 (Mucharz), 2142 kg/m 3 (Radków) and 2448 kg/ m 3 (Tumlin), which is consistent with the density of real sandstones. The modelled samples were compressed by two walls up to failure, both walls moving with a constant velocity 0.2 mm/s. Although this rate is rather higher than that typically used in laboratory uniaxial compression experiments, it is supposed to be sufficiently small to maintain quasi-static conditions in the simulations. The desired speed of the walls was achieved by increasing the velocity gradually (during the initial part of compressing), which also helps to ensure that the sample was loaded quasi-statically without generating an acoustic wave. Full details of the model parameters are presented in Table 2.

Model parameters and comparison between numerical and experimental stressstrain curves
Every specimen represents a medium which, for small external loading, behaves as an ideal elastic body. For larger loads, when some inter-particle bonds break and particles can significantly move away from their initial location, the material exhibits some plasticity. Finally, in large stress concentration regions, particles can separate due to significant stress redistribution when the interaction bonds break, which finally leads to initiation and development of cracks. Visualisation of the numerical samples at the beginning of the loading process and after failure is presented in Table 3, together with a comparison with the laboratory rocks after failure. Basu et al. (2013) carried out profound investigation of failure patterns in different rocks under uniaxial compression. They observed and categorised six possible failure modes, namely, axial splitting, shearing along a single plane, double shear, multiple fracturing, along foliation, and Y-shaped which are schematically represented in Fig. 7. The dominant types of failure after UCS tests in the sandstones tested on the purpose of that article are double shear or Y-shaped depending on the type of rock and strength. Undamaged base surfaces were a characteristic of all samples, which formed cones typical of double shear (Radków and Tumlin) or wedges as in Y-shaped failure (Brenna and Mucharz). The failure modes in numerical tests were quite similar to shearing along a single plane (Brenna, Radków and Tumlin) or double shear (Mucharz).
The relationship between numerical and laboratory samples was found through a process called calibration by comparing numerical and laboratory stress-strain   curves. Parameters of DEM models were calibrated in a way to ensure that the slope of numerical stress-strain curve (Young's modulus E), peak of this curve (UCS) and critical strain e z cr correspond to experimental equivalents. Because the objective of the calibration was to recreate stress-strain curve, only two out of four bond parameters, E b and c b were tuned during the calibration process. Other parameters of bonds like Poisson's ratio v b or internal friction angle φ b were not taken into account due to their trace impact on stress-strain curves (in DEM models); therefore, their values were constant. Such an approach results directly from the algorithm properties, where macroscopic Young's modulus is affected mostly by microscopic Young's modulus and UCS is affected mostly by microscopic cohesion, while the influence of microscopic Poisson's ratio and internal friction angle on the two mentioned macroscopic parameters is negligible or mild. A remarkable advantage of tuning only two parameters was reduction of computational time. Calibration was carried out on the selected laboratory experiments (Brenna -sample BR2, Mucharz -sample MU1, Radków -sample RAD1 and Tumlin -sample TU2). Results are shown in Table 4. Applied calibration method allowed to achieve reasonable consistency between laboratory and simulation values. A comparison between numerical and laboratory stress-strain curves can be found in Fig. 8. Poisson's ratios did not belong to the calibrated parameters; however, its value was also calculated to evaluate deformation of a material in directions perpendicular to the specific direction of loading. Numerically calculated Poisson's ratio of the Brenna and Mucharz were smaller than the same experimental ratios. So, in the case of Brenna, the numerical Poisson's ratio was equal to 0.21 and the experimental ratio was 0.25, and in the case of sandstone Mucharz, it was 0.19 during simulation and 0.25 during the experiment. On the other hand, in the case of sandstones from Radków and Tumlin, higher values were obtained for the numerical samples. Numerically calculated Poisson's ratio of Radków was equal to 0.21 and was almost the same, as the experimental value -0.20. Difference in the case of sandstones from Tumlin was only slightly higher: 0.26 in simulation and 0.22 in experiment. However, it is observed in Table 1 that these values fit well to the value range for the corresponding sandstones and are also very close to the mean values.
Numerically calculated Young's modulus of Brenna was equal to 20.8 GPa and was almost the same as experimental value -20.5 GPa. Similar dependence was obtained for the critical strain: 0.51% for the numerical  The aim of all foregoing efforts was to calibrate the numerical models of sandstones. Numerical and experimental stress-strain curves were overlapped and their characteristic factors were compared. DEM numerical models reproduced selected mechanical parameters similar to the real equivalents. Based on these results, in further studies, calibrated numerical models of sandstones were subjected to different strength trials.

Results
Fracturing of materials such as sandstone is associated with various processes occurring inside the material. In further part of the assay, it was attempted to explain the differences in fracturing process of every sandstone and to investigate what happens inside the sandstone during compression leading to failure. This part was carried out using only DEM simulations. Numerical models used in this investigation were created based on the results derived from the calibration process -it means that the numerical materials (equivalents of the sandstones Brenna, Mucharz, Radków and Tumlin) have the parameters of bonds and particles consistent with Table 2.
Several factors were taken into account to assess the failure process. In the first part, change in the number of bonds was used to follow the appearance of fractures in the material. Then, kinetic energy and its components -kinetic energy of linear motion and kinetic energy of rotational motion -were used to assess the dynamic aspects of the breaking process. Potential energy of bonds was analysed to follow the differences between sandstones in storing energy inside matter. Further, four components of potential energy of bonds were considered (connected with bending, shear, torsion, and tension) to find the dominant interactions between grains of matter. In the final part, all grains were counted into which the sample fell apart.
In the second part, strength of all four numerical sandstones was tested. Real sandstones from Brenna, Mucharz, Radków and Tumlin are commonly used in practical applications; therefore, it means that they are subjected to different external factors threatening the structure of the material. Such factors include breaking forces, various types of impacts, external vibrations and material abrasion. Proposed DEM tests consisted of four experiments, namely, failure resistance test, impact resistance test, vibration resistance test and abrasion resistance test, as a numerical attempt to measure the resistance of sandstones that originated from quarries in Brenna, Mucharz, Radków and Tumlin.

Sandstones under uniaxial compressional regime
Sandstone stress-strain curve can be divided into three main parts: loading phase (before failure), failure and post-failure phase. Within loading stage, it is possible to distinguish the initial phase, elastic phase and plastic deformation. After, dynamic rupture occurs, and then post-failure phase. Characteristics of such a curve provide information on the behaviour of the material during the compression process. As it was shown before, DEM modelling was capable of recreating stress-strain curves of real sandstones (Fig. 8). During numerical loading, different parameters associated with the particles' movements were recorded, and Figs. 9-13 present the compressive process from a microscopic point of view of single particles.
Initially, all particles inside the numerical models were connected by bonds. It means that there were no microcracks or voids in the sense of real material. Together with progressing compression, some bonds started to break apart. In this way, the place with a lacking bond was identical to a fracture appearance. Firstly, the number of bonds (parameter opposite to the number of fractures) was studied as a function of axial displacement (Fig. 9). The number of bonds at every time step was scaled against the initial number of bonds N bonds /N bonds max . Failure of the material was shown as rapid decrease in the number of bonds corresponding to the critical strain of the material e z_cr and stress drop on the stress-strain curve. Sandstones from Radków and Brenna were characterised by similar decrease in the number of bonds (about 40% percent between the initial and final values), whereas this drop was equal to almost 30% in case of numerical Tumlin and only 12% in case of Mucharz. Breaking of numerical Mucharz was also quick and rapid. It seems to be connected with the failure patterns of the samples (Table 3), where numerical Mucharz was the only sample characterised by double shear mode.
To explain this process more profoundly, Fig. 10 presents a comparison between stress-strain curve and number of fractures as a function of axial displacement. For better readability, curves were scaled to maximum values of fractures N frac /N frac max . Each sandstone had similar shape of stress-strain curve, showing a linear increase of stress during the loading stage and rapid failure thereafter. Stress-strain curves were compared against a number of fractures as a function of displacement (dotted line). Because the failure occurred as a single, sudden event, fractures appeared in a similar manner. Numerical sandstones Mucharz and Tumlin were characterised by the appearance of only one single peak in a function describing the number of fractures. Meanwhile, sandstones from Brenna and Radków showed a gentle step before the main failure event (in the range corresponding to the linear part of stress and strain curve). This was a sign of appearance of fractures at the elastic phase.
Particles in the models were connected by bonds, which were, in fact, elastic cylindrical beams, somehow similar to springs (Fig. 3). Due to operating stresses, such beams could store potential energy. Storing and release of potential energy of bonds inside the numerical sandstones under study is presented in the Fig. 11A. Together with progressing compression (loading phase), potential energy was stored in bonds and then sudden release occurred during failure -potential energy was converted into kinetic energy. Peaks of the curves corresponded to peaks of relevant stress-strain curves. All numerical sandstones were characterised by similar shape of the potential energy curves.
Meanwhile, kinetic energy was connected with all possible kinds of movements of particles. In Fig. 11B is presented the dependence between scaled kinetic energy and axial displacement. Kinetic energy behaved in a way opposite to the behaviour of potential energy -it rose significantly after the materials' failure, when potential energy of the bonds was converted into kinetic energy of particle movement and then the kinetic energy decreased. Kinetic energy did not drop to zero because after failure, pieces of broken material were still being compressed. Curves for particular sandstones showed significant similarity, besides the sandstone from Mucharz, which exhibited much more significant decrease in the value of kinetic energy. On comparison with Fig. 9, it is observed that failure in case of sandstone from Mucharz was a quick, rapid event, during which most of the energy was released. For this reason, in  case of Mucharz, kinetic energy drop was much more significant than for other sandstones. Again, it seems to be the result of the fact that numerical Mucharz was characterised by a different failure mode than other samples.
Two types of kinetic energy can be distinguished, linear and rotational, because particles can move from one point to another, as well as they can rotate. The graphs in Fig. 12 give an answer to what was the change of components of kinetic energy: energy of translational movement and energy of rotational movement. The curves again were scaled to maximum value. It is interesting to see that rotations do not play an important role during compressive loading leading to failure, probably even after failure, there were occurring mostly bigger fragments, not single grains able to rotate. Therefore an overwhelming majority belonged to linear kinetic energy, with only very small contribution from rotations.
Similarly, it is possible to divide potential energy of bonds into subtypes, relevant to movements which a cylindrical beam can perform (Fig. 3). Consequently, these components of potential energy of bonds are as follows: potential energy of shear, bending, torsion and tension. From the charts presented in Fig. 13, it can be concluded that for all sandstones the biggest contribution belonged to potential energy of tension/compression and the next one was potential energy of a shearing. Contribution from torsion and bending was very small, almost negligible. This relation was very congenial for every considered numerical sandstone; however, again numerical Mucharz was characterised by slightly different proportion, where the contribution from energy of bending was not so significantly higher in comparison to energy of shearing as in the case of other numerical materials.
After the numerical compression strength test, every numerical sandstone broke into different number of pieces. DEM allowed to calculate precis numbers and weights of the elements (Table 5). Sandstone Brenna fall apart into 19,888 fragments, sandstone Mucharz into 23,384 fragments, sandstone Radków into 33,790 fragments and sandstone Tumlin into 152,953 fragments. Distinguishable case of sandstone Tumlin seemed to originate from the fact that this sandstone consisted of much finer grains than the other sandstones. The initial mass of the samples was: 153 g (Brenna), 113 g (Mucharz), 118 g (Radków) and 135 g (Tumlin). The largest fragments that remained after failure were: 79 g (Brenna), 48 g (Mucharz), 21 g (Radków) and 75 g (Tumlin). Table 5 depicts distribution of the fragments after failure. Mass of the largest fragments was used as a reference to classify all the fragments. Pieces lighter than 10% of the mass of the largest fragment were classified as the smallest grains. Pieces heavier than 90% of the mass of the largest fragment were classified as the heaviest grains. Others were classified as middle grains. An overwhelming majority of the fragments were the smallest ones, mostly weighing less than hundredths and thousandths of a gram. Number of the biggest and middle grains varied from one to three (the biggest) and from four to eight (middle ones). Sandstone Radków stood out from the rest -it had been crushed into a relatively large number of fragments.

Strength of sandstones
Real sandstones used in industrial applications are constantly subjected to various external factors. These factors can be various like the impact by external objects, breakage forces, external vibration or friction, to name a few. In this part, DEM simulations were applied to test how numerical sandstones reacted to external factors during four different tests. For this purpose, rectangular slabs (samples) with dimensions 30 mm × 30 mm of various thicknesses were generated, which were then tested for their strength. The input parameters of the simulation were chosen in accordance with the results obtained in the previous stage of work (during calibration) in an attempt to simulate durability of the numerical sandstones under study -Brenna, Mucharz, Radków and Tumlin.
The proposed tests consisted of four experiments, namely: 1) failure resistance test -an indicator of limit that refers to the power or capacity of a mechanism, structure or system to resist or withstand failure; 2) impact resistance test -measure of the resistance of materials to mechanical impact; 3) vibration resistance test -a measure of the resistance of materials to applied external sinusoidal oscillations and 4) abrasion resistance test -the ability of a surface to resist being worn away by rubbing or friction. Exemplary views of the numerical materials after two different tests are presented in Fig. 14.
The simplified scheme of all four applied tests is presented in Fig. 15. Failure resistance test was carried  out in such a way that the sample on one side was compressed in a vice and the protruding fragment was subjected to external force (Fig. 15A). Material samples with thicknesses 6, 7, 8, 9 and 10 mm were tested. The effect of the test was breaking the protruding piece of the sample. During the test, the applied force was recorded in order to find out the response of each sandstone against the breaking force. The scheme of impact resistance test is presented in the Fig. 15B. Rectangular model of numerical sandstone was tested for its resistance against impact by an external object. The sample was again a 30 mm by 30 mm rectangle with thickness varying from 5 to 10 mm. Above the sample, 5 mm above its upper surface, was located a spherical object with a diameter of 3 mm and density of 11,000 kg/m 3 (close to the density of a lead). The object was given a constant velocity (100 m/s) and made to move towards the sandstone sample. In the third part of testing, the samples were subjected to external vibrations. The scheme of vibration resistance test is presented in Fig. 15C. External oscillations were applied to the same rectangular sample, as during the previous test, with the sample thickness ranging from 5 to 10 mm. Value of the applied oscillations was the same for each sandstone, with the frequency being 8157 Hz (chosen arbitrarily). The last part of the series of resistance tests was the abrasion test (simplified scheme in Fig. 15D). Sample was compressed with a constant force and, at the same time, shearing force was applied. Level of friction was measured by applying bulk friction coefficient factor. The bulk friction coefficient of a sheared material was estimated by measuring the amount of force required to maintain a constant rate of shear at the boundaries. The bulk friction coefficient was defined as the measured shear force divided by the applied normal force.
Results of the tests are presented in Fig. 16. In case of failure resistance test (Fig. 16A) almost linear dependence was noticed between scaled maximum force and thickness of the sample. This dependence was almost the same for the numerical sandstones Mucharz, Radków and Brenna. Results for these three sandstones partially overlapped with each other. Meanwhile, the numerical sandstone Tumlin could withstand much higher breaking force; however, its strength depended much weaker on thickness. It means that the material build of numerical sandstone Tumlin does not gain strength with increasing thickness as significantly as other sandstones.
The effect of the impact resistance experiment is presented in Fig. 16B. As a result of impact, structure of the numerical slab was damaged in the central part of the sample and, consequently, bonds in this area disappeared. Scaled number of fractures as a function of sample thickness was almost constant and independent of thickness. The most numerous cracks appeared for the numerical sandstones Radków and Brenna, almost twice less than for the sandstone Tumlin. In case of numerical sandstone Mucharz, nearly no cracks appeared. It seems that the number of fractures is inversely dependent on UCS (    zone in the centre of the sample was the widest for the sandstone Radków and then, for the sandstone Brenna. Damage zone in case of sandstones Mucharz and Tumlin was much more compact. Damage zone was surrounded by a branching network of cracks, except sandstone Mucharz. In this case, the network of cracks did not appear. The last two tests were vibration resistance test (Fig. 16C) and abrasion resistance test (Fig. 16D). The consequence of vibration resistance test was destruction of the sample structure and its disintegration into pieces. Depending on the type of sandstone and its thickness, the material response to applied external vibrations was different. It is interesting to see that vibration resistance test confirmed the results shown in the Table 5. There was a significant difference between sandstone Radków and other sandstones. Sandstone Radków crumbled into the highest number of crumbs. This value increased with sample thickness; however, this was mild growth. On the opposite, number of grains for other sandstones rather seemed to decrease; additionally, this dependence between the number of grains and sample thickness was almost the same for sandstones Brenna, Mucharz and Tumlin. Meanwhile calculations after abrasion resistance test showed that sandstone Tumlin was characterised by the highest bulk friction coefficient (Fig. 16). Values of this factor for sandstones Brenna, Mucharz and Tumlin were smaller. In this case, bulk friction coefficient seems to directly arise from the size of particles inside the matter. Proposed test showed that sandstone from Tumlin was Figure 16: Results of all four tests: (A) failure resistance test, dependence between scaled maximum force F/F max (force divided by maximum force) and sample thickness; (B) impact resistance test, dependence between scaled number of fractures N frac /N frac max (number of fractures divided by maximum number of fractures) which appeared in the sample after the test and sample thickness; (C) vibration resistance test, number of fragments into which the sample has fallen apart as a function of sample thickness; (D) abrasion resistance test, bulk friction coefficient of the sandstone. characterised by the highest value of friction, together with sandstone from Mucharz. Sandstones from Brenna and Radków, which consisted of relatively bigger grains, were characterised by smaller coefficient of friction. Basu et al. in 2013 reported that rock failure is a serious problem in rock engineering environments and this sentence, despite the passage of years, seems to be still valid. However, recently, computer simulations are becoming a more and more helpful tool in geotechnics together with the increasing power of computers. In the present study numerical method DEM was applied to simulate the behaviour of materials under uniaxial compression. DEM is a technique that is extremely computationally expensive, and the first task of this study was to create the numerical models of dimensions of real samples used in laboratory, which consisted of particles of a size similar to the size of real grains. Four types of sandstones originating from active quarries in Poland (from Brenna, Mucharz, Radków and Tumlin) were used as a pattern. Such sandstones are commonly used in industrial applications. A series of laboratory measurements were carried out on these materials to define their basic properties, basically during uniaxial compression test. Creation of numerical equivalents was a challenging task. Model of a sample with the tiniest grains, Tumlin, consisted of almost 4 million particles. It was necessary to employ the power of the supercomputer to recreate uniaxial compression tests of selected sandstones. In the most problematic case of Tumlin, one test required about 24 hours and 320 processors. Therefore it seems that DEM modelling of real-size samples is still a matter of future.

Summary and Discussion
Macroscopic parameters of real materials were used to calibrate microscopic parameters of numerical models, and calibrated numerical quasi-sandstones were the object of different tests. In the first part of the investigation, it was attempted to explain the differences in the fracturing process for each of the sandstones and to investigate what was happening inside the sandstone during compression leading to failure from the point of view of energies of particles and bonds between them. In the second stage, durability of all four numerical sandstones was tested. Proposed tests consisted of four experiments, namely: failure resistance test -an indicator of limit that refers to the power or capacity of a mechanism, structure or system to resist or withstand failure; impact resistance test -a measure of the resistance of materials to mechanical impact and vibration resistance test -a measure of the resistance of materials to applied external sinusoidal oscillations; abrasion resistance test -the ability of a surface to resist being worn away by rubbing or friction.
However, the proposed approach was full of simplifications, often caused by still limited possibilities of computers. At a material scale, mineralogy and geometric arrangements of particles and voids and/or microcracks together control the mechanical behaviours of rocks. Rock materials contain initial damages, microcracks and so on, and under compression, they undergo a complex process of crack closure and initiation, stable crack growth and unstable cracking and eventually fail (Basu et al. 2013). Meanwhile, the proposed numerical models were calibrated only by comparison of stressstrain curves and calculation of Poisson's ratio. They did not include certain parameters as the initial system of microcracks and/or weakness planes causing strong anisotropy in sedimentary rocks. DEM models should be able to reproduce macroscopic compressive strength, tensile strength (failure patterns as well), dilatancy and scale effect; therefore, an idea for further research can be combining uniaxial compression test with triaxial test to carry out more profound calibration of the models and then introducing to them features responsible for material anisotropy.
At the monoliths' sampling stage, the first field selection was made and the best and most durable monoliths (not weathered, etc.) were selected. Tested samples were cut from these monoliths. The obtained results take into account the scale effect known in the literature, according to which an increase in UCS of samples with small dimensions and slenderness is observed. In situ tests are performed for a larger area (on large-sized samples), where higher loads can be applied, and such tests are more representative of the discontinuities that occur. In general, it is necessary to have very limited confidence in the modelling results and verify the numerical results with laboratory data and the results of measurements of actual structures. In laboratory tests, all parameters are determined on the basis of the adopted models for rock samples. In situ tests also determine all parameters, but for larger samples. However, in numerical analysis, parameters are calculated for a specific construction and changes in the structure under load are monitored.
It seems that at the current stage, numerical modelling of DEM is not able to replace laboratory tests in determining the engineering parameters of materials. Such simulations are still too time consuming, and the calibration is arduous and labour-intensive. In addition, DEM models are not able to effectively recreate real-sized objects. But despite the mentioned limitations of the quantitative approach, DEM works great in the qualitative approach. DEM offers the possibility to control and monitor microscopic parameters and to analyse the breaking processes from the point of view of single particles (their energies, interactions with neighbouring particles or acting forces), leading to better understanding of the phenomenon. An idea for future work is implementation of defects present in a rock matrix (e.g. existing micro cracks, boundaries between crystals) in the model, which strongly affects its mechanical response at least at a laboratory scale. Quantifying defects (type, size, distribution) is very complex to do experimentally, however, possible numerically.