1. bookAHEAD OF PRINT
Journal Details
License
Format
Journal
First Published
09 Nov 2012
Publication timeframe
4 times per year
Languages
English
access type Open Access

Strength of industrial sandstones modelled with the Discrete Element Method

Published Online: 09 Oct 2021
Page range: -
Received: 31 Dec 2020
Accepted: 08 Jul 2021
Journal Details
License
Format
Journal
First Published
09 Nov 2012
Publication timeframe
4 times per year
Languages
English
Abstract

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.

Keywords

Introduction
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 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

The scheme of DEM simulation (after O’Sullivan 2011).

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

Figure 2

The interaction of the two particles by the force–displacement law. The figure is redrawn from Ferdowsi (2014). Particles a and b of radii Ra and Rb are at a distance D with overlap C. Pa and Pb 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.

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

Forces and moments between the particles bonded through rotational elastic-brittle bonds (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: 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 Ks = Ks1 = Ks2 and Kb = Kb1 = Kb2. 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 cb 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 Eb, Poisson's ratio vb, cohesion cb 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 (Fs)max = Fn 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: τij=FτijAij+MtijRijJijσij=-1FσijAij+MbijRijIij\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).

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

Figure 4

Location of the active quarries on the map of Poland from which the tested sandstones originated.

Figure 5

Schematic representation of the uniaxial compression test. Sample is compressed from the top and bottom until failure. On the left – simplified scheme of the test, on the right – laboratory experiment.

Figure 6

Cracked samples after uniaxial compression test.

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 medium-grained 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/L0 where ΔL is the measured displacement and L0 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 ν = εxyz.

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 (Fmax) was determined to the nearest 1 kN. Next, the value of the UCS was calculated as follows: UCS = Fmax/A (MPa) where Fmax is the maximum load at failure (kN) and A is the cross-sectional area (mm2).

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.

Range of results of geomechanical laboratory tests (mean values in brackets).

Lithological typeρ (kg/m3)UCS (MPa)E (GPa)ν (−)
Gogula sandstones from Brenna2391 – 2444 (2429)99 – 107 (102)20.2 – 21.8 (20.8)0.21 – 0.25 (0.22)
Krosno sandstones from Mucharz2627 – 2668 (2655)138 – 153 (143)26.6 – 31.8 (29.2)0.18 – 0.23 (0.19)
Joint sandstones from Radków2119 – 2172 (2142)55 – 60 (58)22.5 – 24.9 (23.4)0.20 – 0.29 (0.24)
Sandstones from Tumlin2355 – 2480 (2448)115 – 128 (123)26.2 – 29.4 (27.4)0.22 – 0.28 (0.25)

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

Detailed parameters of numerical models of sandstones.

Sandstone 1 BrennaSandstone 2 MucharzSandstone 3 RadkówSandstone 4 Tumlin
Sample geometryCylinder:height – 75 mm;radius – 18.75 mmCylinder:height – 75 mm;radius – 18.75 mmCylinder:height – 75 mm;radius – 18.75 mmCylinder:height – 75 mm;radius – 18.75 mm
Particles radii0.25 – 2.0 mm0.25 – 0.5 mm0.25 – 1.0 mm0.125 – 0.5 mm
Particles density2429 kg/m32655 kg/m32142 kg/m32448 kg/m3
Initial number of particles82,266339,324150,6531,155,851
Initial number of bonds297,847821,488500,0243,831,657
Time step1.23×10−5s8.32×10−6s1.06×10−5s5.14×10−6s
Parameters of bonds :- Young's modulus Eb (MPa)- Poisson's ratio vb (dimensionless)- cohesion cb (MPa)- tangent of internal friction angle φb (dimensionless)Eb = 5 284 MPavb = 0.25cb = 13.59 MPatan φb = 1.0Eb= 50 258 MPavb= 0.25cb = 205.33 MPatan φb = 1.0Eb = 12 492 MPavb = 0.25cb = 21.19 MPatan φb = 1.0Eb = 15 144 MPavb = 0.25cb = 44.7 MPatan φb = 1.0
Stroke rate0.2 mm/s0.2 mm/s0.2 mm/s0.2 mm/s
Computational time287 min667 min335 min1440 min
Number of processors9222893320
Model parameters and comparison between numerical and experimental stress–strain 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.

Comparison of failure patterns in DEM modeling (initial model, and model just after the failure - view of bonds between particles) and in laboratory (laboratory failure).

TypeInitial numerical modelNumerical failureLaboratory failure
Brenna sample BR2
Mucharz sample MU1
Radków sample RAD1
Tumlin sample TU2

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

Schematic representation of different possible failure patterns of the sample under uniaxial compression extracted and proposed by Basu et al. (2013).

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, Eb and cb were tuned during the calibration process. Other parameters of bonds like Poisson's ratio vb 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.

Comparison between crucial macroscopic parameters used for model calibration. ‘Sim’ stands for simulation result and ‘lab’ stands for result from the laboratory experiment.

Selected testE (GPa)UCS (MPa)ez cr(%)ν (−)
Brenna sample BR2Sim: 20.8Lab: 20.5Sim: 98Lab: 99Sim: 0.51Lab: 0.52Sim: 0.21Lab: 0.25
Mucharz sample MU1Sim: 29.2Lab: 31.8Sim: 136Lab: 138Sim: 0.47Lab: 0.49Sim: 0.19Lab: 0.23
Radków sample RAD1Sim: 23.4Lab: 22.8Sim: 61Lab: 60Sim: 0.28Lab: 0.29Sim: 0.21Lab: 0.20
Tumlin sample TU2Sim: 27.4Lab: 26.2Sim: 118Lab: 115Sim: 0.46Lab: 0.51Sim: 0.26Lab: 0.22

Figure 8

Comparison between numerical and laboratory stress-strain curves for selected samples: (A) Brenna BR2, (B) Mucharz MU1, (C) Radków RAD1 and (D) Tumlin TU2.

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.

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. 913 present the compressive process from a microscopic point of view of single particles.

Figure 9

Scaled number of bonds Nbonds/Nbonds max (number of bonds divided by maximum number of bonds) inside the sample as a function of axial displacement for all four types of sandstones.

Figure 10

Scaled number of fractures (Nfrac/Nfrac max-number of fractures divided by maximum number of fractures) as a function of axial strain with the stress–strain curve in the background: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.

Figure 11

(A) Scaled potential energy of bonds Epot/Epot max and (B) scaled kinetic energy Ekin/Ekin max for all four types of sandstones.

Figure 12

Comparison between two scaled different components of kinetic energy (Eckin/Eckin max – kinetic energy divided by maximum kinetic energy) plotted as a function of strain: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.

Figure 13

Four different components of scaled potential energy of bonds (Ecpot/Ecpot max – potential energy of bonds divided by maximum potential energy of bonds) plotted together: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.

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 Nbonds/Nbonds 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 Nfrac/Nfrac 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).

Number of grains after failure.

Total number of grainsNumber of the smallest grainsNumber of middle grainsNumber of the biggest grains
Brenna19,88819,87981
Mucharz23,38423,37743
Radków3379033,7641412
Tumlin152,953152,94481

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.

Figure 14

Exemplary presentation of the numerical slab after two different tests: (A) failure resistance test, (B) impact resistance test. View of bonds.

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/m3 (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

Schematic representation of all four tests: (A) failure resistance test, (B) impact resistance test, (C) vibration resistance test and (D) abrasion resistance test. View of the cross-section of the numerical slab.

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

Results of all four tests: (A) failure resistance test, dependence between scaled maximum force F/Fmax (force divided by maximum force) and sample thickness; (B) impact resistance test, dependence between scaled number of fractures Nfrac/Nfrac 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.

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

Network of fractures which appeared in the sample after impact resistance test. View of bonds.

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.

Summary and Discussion

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.

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 stress–strain 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.

Figure 1

The scheme of DEM simulation (after O’Sullivan 2011).
The scheme of DEM simulation (after O’Sullivan 2011).

Figure 2

The interaction of the two particles by the force–displacement law. The figure is redrawn from Ferdowsi (2014). Particles a and b of radii Ra and Rb are at a distance D with overlap C. Pa and Pb 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.
The interaction of the two particles by the force–displacement law. The figure is redrawn from Ferdowsi (2014). Particles a and b of radii Ra and Rb are at a distance D with overlap C. Pa and Pb 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.

Figure 3

Forces and moments between the particles bonded through rotational elastic-brittle bonds (Abe et al. 2014).
Forces and moments between the particles bonded through rotational elastic-brittle bonds (Abe et al. 2014).

Figure 4

Location of the active quarries on the map of Poland from which the tested sandstones originated.
Location of the active quarries on the map of Poland from which the tested sandstones originated.

Figure 5

Schematic representation of the uniaxial compression test. Sample is compressed from the top and bottom until failure. On the left – simplified scheme of the test, on the right – laboratory experiment.
Schematic representation of the uniaxial compression test. Sample is compressed from the top and bottom until failure. On the left – simplified scheme of the test, on the right – laboratory experiment.

Figure 6

Cracked samples after uniaxial compression test.
Cracked samples after uniaxial compression test.

Figure 7

Schematic representation of different possible failure patterns of the sample under uniaxial compression extracted and proposed by Basu et al. (2013).
Schematic representation of different possible failure patterns of the sample under uniaxial compression extracted and proposed by Basu et al. (2013).

Figure 8

Comparison between numerical and laboratory stress-strain curves for selected samples: (A) Brenna BR2, (B) Mucharz MU1, (C) Radków RAD1 and (D) Tumlin TU2.
Comparison between numerical and laboratory stress-strain curves for selected samples: (A) Brenna BR2, (B) Mucharz MU1, (C) Radków RAD1 and (D) Tumlin TU2.

Figure 9

Scaled number of bonds Nbonds/Nbonds max (number of bonds divided by maximum number of bonds) inside the sample as a function of axial displacement for all four types of sandstones.
Scaled number of bonds Nbonds/Nbonds max (number of bonds divided by maximum number of bonds) inside the sample as a function of axial displacement for all four types of sandstones.

Figure 10

Scaled number of fractures (Nfrac/Nfrac max-number of fractures divided by maximum number of fractures) as a function of axial strain with the stress–strain curve in the background: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.
Scaled number of fractures (Nfrac/Nfrac max-number of fractures divided by maximum number of fractures) as a function of axial strain with the stress–strain curve in the background: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.

Figure 11

(A) Scaled potential energy of bonds Epot/Epot max and (B) scaled kinetic energy Ekin/Ekin max for all four types of sandstones.
(A) Scaled potential energy of bonds Epot/Epot max and (B) scaled kinetic energy Ekin/Ekin max for all four types of sandstones.

Figure 12

Comparison between two scaled different components of kinetic energy (Eckin/Eckin max – kinetic energy divided by maximum kinetic energy) plotted as a function of strain: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.
Comparison between two scaled different components of kinetic energy (Eckin/Eckin max – kinetic energy divided by maximum kinetic energy) plotted as a function of strain: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.

Figure 13

Four different components of scaled potential energy of bonds (Ecpot/Ecpot max – potential energy of bonds divided by maximum potential energy of bonds) plotted together: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.
Four different components of scaled potential energy of bonds (Ecpot/Ecpot max – potential energy of bonds divided by maximum potential energy of bonds) plotted together: (A) Brenna, (B) Mucharz, (C) Radków, (D) Tumlin.

Figure 14

Exemplary presentation of the numerical slab after two different tests: (A) failure resistance test, (B) impact resistance test. View of bonds.
Exemplary presentation of the numerical slab after two different tests: (A) failure resistance test, (B) impact resistance test. View of bonds.

Figure 15

Schematic representation of all four tests: (A) failure resistance test, (B) impact resistance test, (C) vibration resistance test and (D) abrasion resistance test. View of the cross-section of the numerical slab.
Schematic representation of all four tests: (A) failure resistance test, (B) impact resistance test, (C) vibration resistance test and (D) abrasion resistance test. View of the cross-section of the numerical slab.

Figure 16

Results of all four tests: (A) failure resistance test, dependence between scaled maximum force F/Fmax (force divided by maximum force) and sample thickness; (B) impact resistance test, dependence between scaled number of fractures Nfrac/Nfrac 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.
Results of all four tests: (A) failure resistance test, dependence between scaled maximum force F/Fmax (force divided by maximum force) and sample thickness; (B) impact resistance test, dependence between scaled number of fractures Nfrac/Nfrac 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.

Figure 17

Network of fractures which appeared in the sample after impact resistance test. View of bonds.
Network of fractures which appeared in the sample after impact resistance test. View of bonds.

Number of grains after failure.

Total number of grainsNumber of the smallest grainsNumber of middle grainsNumber of the biggest grains
Brenna19,88819,87981
Mucharz23,38423,37743
Radków3379033,7641412
Tumlin152,953152,94481

Comparison of failure patterns in DEM modeling (initial model, and model just after the failure - view of bonds between particles) and in laboratory (laboratory failure).

TypeInitial numerical modelNumerical failureLaboratory failure
Brenna sample BR2
Mucharz sample MU1
Radków sample RAD1
Tumlin sample TU2

Range of results of geomechanical laboratory tests (mean values in brackets).

Lithological typeρ (kg/m3)UCS (MPa)E (GPa)ν (−)
Gogula sandstones from Brenna2391 – 2444 (2429)99 – 107 (102)20.2 – 21.8 (20.8)0.21 – 0.25 (0.22)
Krosno sandstones from Mucharz2627 – 2668 (2655)138 – 153 (143)26.6 – 31.8 (29.2)0.18 – 0.23 (0.19)
Joint sandstones from Radków2119 – 2172 (2142)55 – 60 (58)22.5 – 24.9 (23.4)0.20 – 0.29 (0.24)
Sandstones from Tumlin2355 – 2480 (2448)115 – 128 (123)26.2 – 29.4 (27.4)0.22 – 0.28 (0.25)

Detailed parameters of numerical models of sandstones.

Sandstone 1 BrennaSandstone 2 MucharzSandstone 3 RadkówSandstone 4 Tumlin
Sample geometryCylinder:height – 75 mm;radius – 18.75 mmCylinder:height – 75 mm;radius – 18.75 mmCylinder:height – 75 mm;radius – 18.75 mmCylinder:height – 75 mm;radius – 18.75 mm
Particles radii0.25 – 2.0 mm0.25 – 0.5 mm0.25 – 1.0 mm0.125 – 0.5 mm
Particles density2429 kg/m32655 kg/m32142 kg/m32448 kg/m3
Initial number of particles82,266339,324150,6531,155,851
Initial number of bonds297,847821,488500,0243,831,657
Time step1.23×10−5s8.32×10−6s1.06×10−5s5.14×10−6s
Parameters of bonds :- Young's modulus Eb (MPa)- Poisson's ratio vb (dimensionless)- cohesion cb (MPa)- tangent of internal friction angle φb (dimensionless)Eb = 5 284 MPavb = 0.25cb = 13.59 MPatan φb = 1.0Eb= 50 258 MPavb= 0.25cb = 205.33 MPatan φb = 1.0Eb = 12 492 MPavb = 0.25cb = 21.19 MPatan φb = 1.0Eb = 15 144 MPavb = 0.25cb = 44.7 MPatan φb = 1.0
Stroke rate0.2 mm/s0.2 mm/s0.2 mm/s0.2 mm/s
Computational time287 min667 min335 min1440 min
Number of processors9222893320

Comparison between crucial macroscopic parameters used for model calibration. ‘Sim’ stands for simulation result and ‘lab’ stands for result from the laboratory experiment.

Selected testE (GPa)UCS (MPa)ez cr(%)ν (−)
Brenna sample BR2Sim: 20.8Lab: 20.5Sim: 98Lab: 99Sim: 0.51Lab: 0.52Sim: 0.21Lab: 0.25
Mucharz sample MU1Sim: 29.2Lab: 31.8Sim: 136Lab: 138Sim: 0.47Lab: 0.49Sim: 0.19Lab: 0.23
Radków sample RAD1Sim: 23.4Lab: 22.8Sim: 61Lab: 60Sim: 0.28Lab: 0.29Sim: 0.21Lab: 0.20
Tumlin sample TU2Sim: 27.4Lab: 26.2Sim: 118Lab: 115Sim: 0.46Lab: 0.51Sim: 0.26Lab: 0.22

Abe, S., V. Boros, W. Hancock, and D. Weatherley (2014), ESyS-particle tutorial and user's guide. Version 2.3.1, https://launchpad.net/esys-particle.AbeS.BorosV.HancockW.WeatherleyD.2014ESyS-particle tutorial and user's guideVersion 2.3.1, https://launchpad.net/esys-particle.Search in Google Scholar

Basu, A., D. A. Mishra, and K. Roychowdhury (2013), Rock failure modes under uniaxial compression, Brazilian, and point load tests, Bull Eng Geol Environ, 72:457–475, DOI 10.1007/s10064-013-0505-4BasuA.MishraD. A.RoychowdhuryK.2013Rock failure modes under uniaxial compression, Brazilian, and point load testsBull Eng Geol Environ7245747510.1007/s10064-013-0505-4Open DOISearch in Google Scholar

Carmona, H.A., F.K. Wittel, F. Kun, and H.J. Herrmann (2008), Fragmentation processes in impact of spheres, Phys. Rev. E 77, 5, 051302, DOI: 10.1103/PhysRevE.77.051302.CarmonaH.A.WittelF.K.KunF.HerrmannH.J.2008Fragmentation processes in impact of spheresPhys. Rev. E77505130210.1103/PhysRevE.77.051302Open DOISearch in Google Scholar

Courtney, T.H. (1990), Mechanical Behavior of Materials, McGraw-Hill, New York.CourtneyT.H.1990Mechanical Behavior of MaterialsMcGraw-HillNew YorkSearch in Google Scholar

Cundall, P.A. (1971), A computer model for simulating progressive, large-scale movement in blocky rock systems, Proc. Symp. Int. Soc. Rock Mech. 2, 8.CundallP.A.1971A computer model for simulating progressive, large-scale movement in blocky rock systemsProc. Symp. Int. Soc. Rock Mech.28Search in Google Scholar

Cundall, P.A. (1974), A computer model for rock-mass behavior using interactive graphics for the input and output of geometrical data, National Technical Information Service, Report no. AD/A-001 602.CundallP.A.1974A computer model for rock-mass behavior using interactive graphics for the input and output of geometrical data, National Technical Information Service, Report no. AD/A-001 602Search in Google Scholar

Cundall, P.A., and O.D.L. Strack (1979), A discrete numerical model for granular assemblies, Geotechnique 29, 1, 47–65, DOI: 10.1680/geot.1979.29.1.47.CundallP.A.StrackO.D.L.1979A discrete numerical model for granular assembliesGeotechnique291476510.1680/geot.1979.29.1.47Open DOISearch in Google Scholar

Do, H. Q., Aragón, A. M., and Schott, D. L. (2017). Automated discrete element method calibration using genetic and optimization algorithms. In EPJ Web of Conferences (Vol. 140). EDP Sciences. https://doi.org/10.1051/epjconf/201714015011.DoH. Q.AragónA. M.SchottD. L.2017Automated discrete element method calibration using genetic and optimization algorithmsInEPJ Web of Conferences140EDP Sciences. https://doi.org/10.1051/epjconf/201714015011.Search in Google Scholar

Egholm, D.L. (2007), A new strategy for discrete element numerical models: 1. Theory, J. Geophys. Res. 112, B5, B05203, DOI: 10.1029/2006JB004557.EgholmD.L.2007A new strategy for discrete element numerical models: 1Theory, J. Geophys. Res.112B5B0520310.1029/2006JB004557Open DOISearch in Google Scholar

Ferdowsi, B. (2014), Discrete element modeling of triggered slip in faults with granular gouge: application to dynamic earthquake triggering, doctoral dissertationFerdowsiB.2014Discrete element modeling of triggered slip in faults with granular gouge: application to dynamic earthquake triggering, doctoral dissertationSearch in Google Scholar

Fraige, F.Y., and P.A. Langston (2004), Integration schemes and damping algorithms in distinct element models, Adv. Powder Technol. 15, 2, 227–245, DOI: 10.1163/156855204773644454.FraigeF.Y.LangstonP.A.2004Integration schemes and damping algorithms in distinct element modelsAdv. Powder Technol.15222724510.1163/156855204773644454Open DOISearch in Google Scholar

Griffith, A.A. (1921), The phenomena of rupture and flow in solids, Philosophic. Trans. Roy. Soc. London A, 221, 582–593, DOI: 10.1098/rsta.1921.0006.GriffithA.A.1921The phenomena of rupture and flow in solidsPhilosophic. Trans. Roy. Soc. London A22158259310.1098/rsta.1921.0006Open DOISearch in Google Scholar

Hertzberg, R.W. (1976), Deformation and Fracture Mechanics of Engineering Materials, Wiley, New York.HertzbergR.W.1976Deformation and Fracture Mechanics of Engineering MaterialsWileyNew YorkSearch in Google Scholar

Jaeger, J.C., N.G.W. Cook, and R. Zimmerman (2007), Fundamentals of Rock Mechanics, 4th ed., Blackwell Publ., Malden, 488 pp.JaegerJ.C.CookN.G.W.ZimmermanR.2007Fundamentals of Rock Mechanics4th ed.Blackwell Publ.Malden488Search in Google Scholar

Kazerani, T. (2013), A discontinuum-based model to simulate compressive and tensile failure in sedimentary rock, J. Rock Mech. Geotech. Eng. 5, 5, 378–388, DOI: 10.1016/j.jrmge.2013.07.002.KazeraniT.2013A discontinuum-based model to simulate compressive and tensile failure in sedimentary rockJ. Rock Mech. Geotech. Eng.5537838810.1016/j.jrmge.2013.07.002Open DOISearch in Google Scholar

Klejment P. (2020), The Microscopic Insight into Fracturing of Brittle Materials with the Discrete Element Method (Doctoral Dissertation), Publications of the Institute of Geophysics Polish Academy of Sciences – Geophysical Data Bases, Processing and Instrumentation 427 (A-31), DOI: 10.25171/InstGeoph_PAS_Publs-2020-001KlejmentP.2020The Microscopic Insight into Fracturing of Brittle Materials with the Discrete Element Method (Doctoral Dissertation)Publications of the Institute of Geophysics Polish Academy of Sciences – Geophysical Data Bases, Processing and Instrumentation 427 (A-31),10.25171/InstGeoph_PAS_Publs-2020-001Open DOISearch in Google Scholar

Klemm A., D. Wiggins (2016), Sustainability of Construction Materials. Second edition, Woodhead Publishing, Cambridge.KlemmA.WigginsD.2016Sustainability of Construction MaterialsSecond editionWoodhead PublishingCambridgeSearch in Google Scholar

Li, Y.-G., ed. (2012), Imaging, Modeling and Assimilation in Seismology, De Gruyter, Berlin.LiY.-G.ed.2012Imaging, Modeling and Assimilation in SeismologyDe GruyterBerlinSearch in Google Scholar

Łukaszewski, P. (2003), Development of fracture processes in Silesian Carboniferous sandstones, Geological Quarterly, vol. 47, (1), 29–38.ŁukaszewskiP.2003Development of fracture processes in Silesian Carboniferous sandstonesGeological Quarterly4712938Search in Google Scholar

Łukaszewski, P. (2013), The deformation of flysch sandstones in a complex state of stress, Wydawnictwo Uniwersytetu Warszawskiego, 221, (in Polish, with English summary).ŁukaszewskiP.2013The deformation of flysch sandstones in a complex state of stressWydawnictwo Uniwersytetu Warszawskiego221(in Polish, with English summary).Search in Google Scholar

Mora, P., and D. Place (1993), A lattice solid model for the nonlinear dynamics of earthquakes, Int. J. Modern Phys. C 4, 6, 1059–1074, DOI: 10.1142/S0129183193000823.MoraP.PlaceD.1993A lattice solid model for the nonlinear dynamics of earthquakesInt. J. Modern Phys. C461059107410.1142/S0129183193000823Open DOISearch in Google Scholar

Mora, P., and D. Place (1994), Simulation of the frictional stick-slip instability, Pure Appl. Geophys. 143, 61–87, DOI: 10.1007/BF00874324.MoraP.PlaceD.1994Simulation of the frictional stick-slip instabilityPure Appl. Geophys.143618710.1007/BF00874324Open DOISearch in Google Scholar

Mora, P. and D. Place (2002). “Stress correlation function evolution in lattice solid elasto-dynamic models of shear and fracture zones and earthquake prediction.” Pure and Applied Geophysics, 159(10): 2413–2427.MoraP.PlaceD.2002“Stress correlation function evolution in lattice solid elasto-dynamic models of shear and fracture zones and earthquake prediction.”Pure and Applied Geophysics1591024132427Search in Google Scholar

Munjiza, A. (2004), The Combined Finite-Discrete Element Method, John Wiley & Sons, Chichester.MunjizaA.2004The Combined Finite-Discrete Element MethodJohn Wiley & SonsChichesterSearch in Google Scholar

Nitka, M., and J. Tejchman (2015), Modelling of concrete behaviour in uniaxial compression and tension with DEM, Granular Matter 17, 1, 145–164, DOI: 10.1007/s10035-015-0546-4.NitkaM.TejchmanJ.2015Modelling of concrete behaviour in uniaxial compression and tension with DEMGranular Matter17114516410.1007/s10035-015-0546-4Open DOISearch in Google Scholar

O’Sullivan, C. (2011), Particulate Discrete Element Modelling: A Geomechanics Perspective, Spon Press/Taylor and Francis, London.O’SullivanC.2011Particulate Discrete Element Modelling: A Geomechanics PerspectiveSpon Press/Taylor and FrancisLondonSearch in Google Scholar

O’Sullivan, C., and J.D. Bray (2004), Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme, Eng. Comput. 21, 2–4, 278–303, DOI: 10.1108/02644400410519794.O’SullivanC.BrayJ.D.2004Selecting a suitable time step for discrete element simulations that use the central difference time integration schemeEng. Comput.212–427830310.1108/02644400410519794Open DOISearch in Google Scholar

Onate, E., and J. Rojek (2004), Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems, Computer Meth. Appl. Mech. Eng. 193, 27–29, 3087–3128, DOI: 10.1016/j.cma.2003.12.056.OnateE.RojekJ.2004Combination of discrete element and finite element methods for dynamic analysis of geomechanics problemsComputer Meth. Appl. Mech. Eng.19327–293087312810.1016/j.cma.2003.12.056Open DOISearch in Google Scholar

Orowan, E. (1949), Fracture and strength of solids, Rep. Prog. Phys. 12, 1, 185–232, DOI: 10.1088/0034-4885/12/1/309.OrowanE.1949Fracture and strength of solidsRep. Prog. Phys.12118523210.1088/0034-4885/12/1/309Open DOISearch in Google Scholar

Potyondy, D.O., and P.A. Cundall (2004), A bonded-particle model for rock, Int. J. Rock Mech. Min. Sci. 41, 8, 1329–1364, DOI: 10.1016/j.ijrmms.2004.09.011.PotyondyD.O.CundallP.A.2004A bonded-particle model for rockInt. J. Rock Mech. Min. Sci.4181329136410.1016/j.ijrmms.2004.09.011Open DOISearch in Google Scholar

Rojek, J. (2007), Modelowanie i symulacja komputerowa złożonych zagadnień mechaniki nieliniowej metodami elementów skończonych i dyskretnych, Prace IPPT PAN, Warszawa (in Polish).RojekJ.2007Modelowanie i symulacja komputerowa złożonych zagadnień mechaniki nieliniowej metodami elementów skończonych i dyskretnychPrace IPPT PANWarszawa (in Polish)Search in Google Scholar

Sochor, M. (1998), Strength of Materials I, CTU Publishing House, Prague.SochorM.1998Strength of Materials ICTU Publishing HousePragueSearch in Google Scholar

Wang Y., S. Abe, S. Latham and, P. Mora (2006), Implementation of Particle-scale Rotation in the 3-D Lattice Solid Model, Pure and Applied Geophysics 163, 1769–1785, doi 10.1007/s00024-006-0096-0WangY.AbeS.LathamS.MoraP.2006Implementation of Particle-scale Rotation in the 3-D Lattice Solid ModelPure and Applied Geophysics1631769178510.1007/s00024-006-0096-0Open DOISearch in Google Scholar

Wang, Y.C. (2009), A new algorithm to model the dynamics of 3-D bonded rigid bodies with rotations, Acta Geotech. 4, 2, 117–127, DOI: 10.1007/s11440-008-0072-1.WangY.C.2009A new algorithm to model the dynamics of 3-D bonded rigid bodies with rotationsActa Geotech.4211712710.1007/s11440-008-0072-1Open DOISearch in Google Scholar

Wang, Y.C., S. Xue, and J. Xie (2012), Discrete element method and its applications in earthquake and rock fracture modelling. In: Y.-G. Li (ed.) Imaging, Modeling and Assimilation in Seismology, De Gruyter, 235–262.WangY.C.XueS.XieJ.2012Discrete element method and its applications in earthquake and rock fracture modellingIn:LiY.-G.(ed.)Imaging, Modeling and Assimilation in SeismologyDe Gruyter235262Search in Google Scholar

Weatherley, D.K., V.E. Boros, W.R. Hancock, and S. Abe (2010), Scaling benchmark of EsyS-Particle for elastic wave propagation simulations, 2010 IEEE Sixth Int. Conf. e-Science, Brisbane, Australia, 277–283, DOI: 10.1109/eScience.2010.40.WeatherleyD.K.BorosV.E.HancockW.R.AbeS.2010Scaling benchmark of EsyS-Particle for elastic wave propagation simulations2010 IEEE Sixth Int. Conf. e-ScienceBrisbane, Australia27728310.1109/eScience.2010.40Open DOISearch in Google Scholar

Recommended articles from Trend MD

Plan your remote conference with Sciendo