Strength of industrial sandstones modelled with the Discrete Element Method

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 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 1971, 1974, 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 1971, 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 large-scale 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.

Figure 1

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, 1994). The particle rotation is governed by the Euler's equations.

Figure 2

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 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).

Figure 3

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: (1) Fr=Kr ΔrFs1=Ks1 Δs1Fs2=Ks2 Δs2Mt=Kt ΔαtMb1=Kb1 Δαb1Mb2=Kb2 Δαb2 \matrix{{{F_r}} & = & {{K_r}\,\Delta r} \cr {{F_{s1}}} & = & {{K_{s1}}\,\Delta {s_1}} \cr {{F_{s2}}} & = & {{K_{s2}}\,\Delta {s_2}} \cr {{M_t}} & = & {{K_t}\,\Delta {\alpha _t}} \cr {{M_{b1}}} & = & {{K_{b1}}\,\Delta {\alpha _{b1}}} \cr {{M_{b2}}} & = & {{K_{b2}}\,\Delta {\alpha _{b2}}}} where F represents forces, M torques, K rigidity constant, Δr and Δs are relative translational displacement and Δα 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 b - bending 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: (2) τij=FτijAij+Mtij⋅RijJijσij=-1⋅FσijAij+Mbij⋅RijIij \matrix{{{\tau _{ij}} = {{{F_{\tau ij}}} \over {{A_{ij}}}} + {M_{tij}} \cdot {{{R_{ij}}} \over {{J_{ij}}}}} \cr {{\sigma _{ij}} = - 1 \cdot {{{F_{\sigma ij}}} \over {{A_{ij}}}} + {M_{bij}} \cdot {{{R_{ij}}} \over {{I_{ij}}}}}} 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 non-linear 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 (1971, 1974) and Cundall and Strack (1979) and further explanations can be found in Rojek (2007). Details of DEM implementation in ESyS-Particle are given in Wang et al. (2006, 2012), Wang (2009), Abe et al. (2014, Ferdowsi (2014) and Weatherley et al. (2010).

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₀ where ΔL is the measured displacement and L₀ 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 length-to-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: UCS = F_max /A (MPa) where F_max is the maximum load at failure (kN) and A is the cross-sectional area (mm²).

During the tests, experimental deformation stress–strain 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. The highest values of strength (138–153 MPa) were obtained for medium-grained Krosno sandstone from Mucharz, while the lowest (55–60 MPa) values were obtained for coarse-grained joint sandstone from Radków. In the case of coarse- and medium-grained Godula sandstone from Brenna and medium-grained sandstone from Tumlin, intermediate UCS values (99–107 and 115–128 MPa) were obtained. The values of the Young's modulus (E) were also varied. For the examined sandstone samples – under uniaxial compression conditions – the Young's modulus E varied from 20.2 to 31.8 GPa (Table 1). The lowest values (20.2–21.8 GPa) were characteristic for coarse- and medium-grained Godula sandstone from Brenna. The highest values, 26.2–29.4 and 26.6–31.8 GPa, 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.

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³ (Brenna), 2655 kg/m³ (Mucharz), 2142 kg/m³ (Radków) and 2448 kg/m³ (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.

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).

Figure 7

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 ε_{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.

Figure 8

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 model and 0.52% in the experiment. Difference in the case of UCS was only slightly higher: 98 MPa in simulation and 99 MPa in experiment. Young's modulus of the numerical sample of Mucharz was equal to 29.2 GPa and was close to the experimental value – 31.8 GPa. Critical strain was also very similar: 0.47% in simulation and 0.49% in experiment. Difference in case of UCS was equal to 136 MPa during simulation and 138 MPa during the experiment. Numerical Young's modulus of Radków was equal to 23.4 GPa, with the experimental value being 22.8 GPa. Critical strain was also very similar: 0.28% in simulation and 0.29% in experiment. Difference in the case of UCS was 61 MPa during simulation and 60 MPa during the experiment. Numerical Young's modulus of Tumlin was equal to 27.4 GPa and the experimental value was equal to 26.2 GPa. However, in case of critical strain, the difference was more significant: 0.46% in simulation and 0.51% in experiment. Difference in the case of UCS was as follows: 118 MPa during simulation and 115 MPa during the experiment.

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.

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.

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

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 ε_{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.

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.

Figure 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³ (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.

Figure 15

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.

Figure 16

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 (Table 4). Sandstone with the highest UCS (Mucharz – 136 MPa) was characterised by the smallest number of fractures. Also, Tumlin – 118 MPa, Brenna – 98 MPa, and Radków – 61 MPa had the number of fractures increasing with the lowering value of UCS. Visualisation of the networks of cracks after the impact is presented in Fig. 17 which shows a top view of the samples after the test for the boundary values of thickness. Depending on the material type, different network of cracks appeared.

Figure 17

Each type of sandstone had its own specific pattern according to which a network of cracks appeared. Damage 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 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.