## 1. Introduction

The present chapter describes numerical simulation of sand geomaterial with the discrete element method (DEM), based on input data acquired from physical experiments. First of all, sandy soil is composed of different shape and size particles. That means, it is necessary to provide sieving test to determine granulometry curve of investigated sand. After this step, it is used scanning electronic microscope (SEM) to determine main sand particles morphology parameters, namely area, perimeter, form coefficient, angularity, roundness, circularity, and sphericity coefficients, respectively. Next physical experiments level—sand compression with oedometer. Having all this information, it is possible to start numerical DEM modelling of sand compression test.

Granulometric curve and separate particle shape discretization are one of the most relevant existing problems for numerical modelling via DEM. There are no direct recommendations for a single‐particle subscription with spheres (single‐particle subscription level based on spheres size and quantity). Subscribing soil particles with DEM, it is very important to choose size, shape, and physical properties of modelled particles. It is possible to model realistic size oedometer device without not scaled particles if quantity of modelled particles is small. Modelling of soil which is subscribed with a big quantity of particles is not possible because of the limitations of computer calculation capacity. All simplifications are accepted in order to decrees the calculation time. Another relevant problem, arising from all accepted simplifications, is a very small quantity of the authors working on experimental and numerical testing validation. This fact only proves that there are still a lot of problems related to dispersive systems, which are modelled with DEM.

## 2. Sand composition and morphology parameters

The investigated area is located in the southern part of Lithuanian continental coastal zone of Baltic Sea, the northern part of Klaipeda city [1]. To the north of the Klaipeda city, only the immediate near‐shore contains a sandy strip of Holocene marine sediments (m IV), which occurs in up to 4–5 m of the sea water depth. The material composing the near‐shore sediments in the continental coast zone mainly consists of different sand, where prevailing medium coarse and fine sand [2, 3] with admixture of gravel and organic matter [4, 5]. This sand was used for investigations (**Figure 1**). The coordinates of sampling location [1] are 55°46’4.07″, 21°4’39.06″ (WGS).

The average density of particles (*ρ _{s}*) value of marine sands is 2.67 Mg/m

^{3}and varies from 2.65 to 2.71 Mg/m

^{3}, respectively. The bulk density of the sand varies from 1.83 to 2.09 Mg/m

^{3}, where average is 1.98 Mg/m

^{3}. Regardless to the genesis of marine sand and their grain size distribution, the mineral composition consists mainly of quartz. The natural moisture content depends on the degree of saturation of water and ranges from 13.7 to 27.7%. The void ratio (

*e*) in fine sand varies from 0.474 to 0.778, respectively [4–6]. The recent marine sediments (m IV) were formed in the coastal zone; therefore, a distinctive morphological feature of grain shape has high sphericity, where

*P*= 0.84 [6]. Investigated sand mineralogical composition determined in reference [7], where basic consistency of the sand is ∼85% quarts, ∼6% feldspar with remaining contribution of carbonate, mica and some other minerals. The grading curve of investigated soil is given in

**Figure 2**.

Thirty‐three different sand grains corresponding nine sand fractions [8] have been examined, the total number of examined grains being 297. Investigation of grains morphology parameters provided with SEM and view analysis program. Example of investigation process is given in **Figure 3**. Smaller fractions within the 0.0063–0.15 mm are omitted in the panoramic view due to the lack of space in SEM camera. Rest of all fractions was investigated separately.

The main morphological parameters of investigated sand fractions employed in 2D view analysis [9, 10] include these parameters, namely area (mm^{2}), equivalent diameter (mm), sphericity, circularity, form coefficient, angularity. Final decision to investigate particles only in 2D was accepted according to literature analysis.

2D and 3D view analysis and determination of grain sphericity results accuracy 5–10% [11, 12]. Due to this fact, it is saved a lot of time which is spent for morphological parameters investigations with SEM and view analysis program. Small results differences between 2D and 3D analysis are due to particle landing on the investigation table with the biggest particle gravity centre (**Figure 4**). Particle stays on investigation table with the biggest stability surface in contact with investigation table [13]. In this case, particles maximum length is in a parallel line with investigation table surface and particles height which is perpendicular for investigation tables is the smallest. Results inaccuracy can be increased if it is investigated mainly flat grains [14].

The analysis of the determined morphological parameters within each fraction shows that the shape of all grains is sufficiently similar—It tests the grains differ principally only in size. The obtained mean 2D morphological parameters of the particles are given in the **Table 1**.

Morphological parameter | Mean value |
---|---|

Area (mm^{2}) | 0.112 |

Equivalent diameter (mm) | 0.340 |

Sphericity | 0.836 |

Circularity | 0.515 |

Form coefficient | 0.702 |

Angularity | 0.410 |

Particle shape |

The mean shape of investigated sand grain was determined using [16, 17] given solutions for sand shape characterization according to the particle sphericity and roundness. In this case, it is not necessary to use Fourier descriptors [18]. The analysis of the change of morphology parameters of investigated Baltic Sea soil‐type according to the equivalent diameter increment is shown in **Figure 5**.

## 3. Experimental sand compression with oedometer

Compression tests were performed using computer controlled (fully automatic) oedometer test device. The same sand soil used for morphology parameters investigation was examined in compression tests with the standard oedometer device (ring height *H* = 3.39 cm; ring diameter *D* = 7.14 cm). Testing program consists of the study of the influence of vertical stress ramp (**Figure 6**) on the compression results, the study of crashing of separate particles, and upgrade of testing procedure, respectively.

Obtained results, which are presented in **Figure 6**, clearly show that loading ramp is not having any effect for soil deformation. Due to this reason, it was accepted to use 400 kPa/min loading ramp for rest of all tests with oedometer. From this point, it is necessary to check two things before numerical oedometer modelling:

Does the examined particles having any crushing effect with used maximum vertical stress 400 kPa on provided oedometer tests.

Where are the errors on computer controlled oedometer testing procedure, because in

**Figure 6**given results start not from 0 kPa.

The solution for the checking of particles crushing is to provide a series of oedometer tests with different maximum vertical stress. After each test, it is necessary to run sieving test and to check granulometry change. If the granulometry composition is not changing, that means, particles with accepted maximum vertical stress loading are not crushed. Small change of mass after sieving test can be revealed, but this happens due to particles loss during sieving procedure. Normally, particle loss can be up to 1%. Obtained sieving tests results after different maximum loading are given in **Table 2**.

Knowing that tested soil particles cannot have crushing effect (the analysis of obtained results revealed that sand particles are not crushed up to 800 kPa maximum vertical loading) makes more simple DEM oedometer test modelling.

The main finding from the compression test procedure is presented in **Figure 7**. Comparison of the standard and improved procedures consists of two different stages. The first comparison stage involves vertical loading values from 0 to 100 kPa. In this stage, the standard test procedure cannot evaluate immediate settlements [20], since for contact ensurance between porous stone and soil, 10 kPa vertical stress is established. Therefore, **Figure 7** shows results not from 0 kPa.

The improved testing procedure is better than the standard one, as porous stone is positioned in the calibrated height, selected according to the initial sample height. Performing compression test with the improved testing procedure, it is possible to analyse compression curve from 0 kPa.

The standard and improved testing procedures differ in their obtained compression curves values only, and the character of their vertical strain values remains the same. In the standard procedure, a smaller compressed sample deformation (*ε* = 1.6168%) is obtained comparing with the improved testing procedure (*ε* = 3.1763%). The comparison of different compression procedures allow the judicious selection of DEM compression modelling as the preferred methodology since it is an improved testing methodology and allows the analysis of compression results from 0 kPa.

The soil compression curve mostly obtained in literature is given in **Figure 8** [21].

Analysing soil compression curve given in **Figure 8**, it comes obvious that soil density or void ratio for low stress values does not change. Void ratio change is obtained only when soil is loaded with high stress values. Such interpretation of the compression curve is not reliable due to two reasons:

In the contact place of porous stone and soil sample, contact settlements occur at the low stress due to not perfect soil surface. Contact settlements ensure good porous stone and soil sample contact between each other.

When testing clay soils with large sand or gravel particles, sample surface imperfections occur during soil sample preparation on the sample side, in a contact with the oedometer ring surface. Due to these surface imperfections, the soil sample expands in the horizontal direction during the loading process, and it continues until a sufficient contact with the oedometer ring. This can be considered as an explanation of the appearance of additional sample settlement.

If the soil sample has loading, unloading and reloading steps, **Figure 8** given compression curve is suitable only for the reloading step. Other authors present soil compression curves with vertical loading from 5 to 10 kPa [22–26]. In this case, it is not necessary to show what happens with the soil compression curve when vertical stress is from 0 to 5 or 10 kPa.

## 4. Numerical simulation of compression test with DEM

Evaluation of the actual compressibility properties via soil compression tests is important for employment of subsequent numerical analysis of stress and strain state of ground subjected by supplement loading (e.g. loads transmitted via foundations from a superstructure, interaction of structure and soil strata, etc.). The confined compression (oedometer) test is approved as a relatively fast and simple laboratory test. It is performed under different conditions, loading paths and durability. Test conditions depend on physical changes in multiphase system of soil, generally related with reorganisation of soil grains, that of initial change of skeleton in cohesive soils and velocity of water filtration for saturated soils. One must note that in some cases the duration of testing procedure for prediction of long‐term soil behaviour and in other specific cases is very long [25], thus sometimes taking into account the time and test cost ratio it is not worth even to start the test. On the other hand, one cannot qualitatively explain the variation of soil compression test results, basing on some processing of already known testing and analysis data separately or in concern with view analysis. It is obvious that having not identified the actual physical mechanism for soil grain reorganisation during compression process and its peculiarities, one cannot explain the observed scatter of results under disposal. This mechanism cannot be recorded by applying usual techniques of testing and view analysis but can be simulated applying the relevant DEM techniques (initiated by [27–29]): that of the mathematical models of processes and discrete models for soil grains [1].

The above‐mentioned circumstances, as well as the permanently reducing computational costs, the development of numerical techniques and software in the field of multi‐scale analysis (including the particle strata mechanics), initiated a fast development and the applications numerical analysis in the field of the soil behaviour by means of the DEM. Such approach, combined with an experimental analysis for validation and calibration of the mathematical models, is definitely a promising one, allowing to reduce the price and quantity of laboratory tests and a reasonless conservatism in determining the values of mechanical properties of soils in near future.

Numerical simulation via DEM was provided using DEMMAT code [30]. This original program called DEMMAT has been written in Compaq Visual FORTRAN language. The quality of implementation is handled by a physically observable behaviour of interactions: particle–particle, particle–wall, particle–bottom and/or top plate (porous stone), and by the validation with the results obtained from physical experiments. The application of DEM to a system involves the following basic steps [30]:

Set‐u of initial conditions of the particles and the walls, top and bottom plate (porous stone);

Searching for the contacts between particles;

Applying interaction laws (calculating forces and moments) to all particle–particle, particle–wall, particle–bottom and/or top plate (porous stone);

Applying Newton's second law to determine particle motion. Predicting positions, velocities, accelerations and higher‐order time derivatives;

Updating the current state of particles and the walls;

Assigning a new time

*t = t + Δt*;If current time is within entire time period (

*t ≤ T*)—go to the step 2, otherwise go to the next step;Post‐processing and visualization.

More detailed information about used models and parameters in DEM simulation is in references [1, 30].

Numerical DEM simulations are performed with the same granulometric curve as obtained in experimental testing for the Baltic Sea sand. The used numerical granulometric curve with mean experimentally determined particle shape is given in **Figure 9**.

In this study, recalculation of time period is based on reference [31]:

where *E _{p}*—modelled particle deformation modulus, MPa;

*v*—loading velocity, m/s

^{2};

*v*—initial particle velocity, m/s;

_{ij}*ρ*—particle density, kg/m

^{3};

*R*—smallest modelled particle radius, m.

_{min}Seeking to increase the calculation accuracy, it is necessary to re‐calculate the time period, thereby obtaining the real‐time period for computation:

Particles position and contact forces change according to *Δt* are given in **Figure 10**.

The test process and parameters for simulations are presented in **Table 3**. The actual Young's modulus *E _{oed}* = 200 GPa for oedometer volume parts is employed for simulations with Poisson’s ratio

*υ*= 0.3 and density

*ρ*= 7850 kg/m

^{3}, and simulated oedometer height

*H*= 0.00484 m and diameter

*D*= 0.0102 m. During simulations reached maximum strain was 1.76% (as in the experimental testing). Accepted simulated particles Poisson’s ratio

*υ*= 0.17, oedometer walls Poisson’s ratio

*υ*= 0.3.

Quantity | Symbol | Unit | Numerical simulation | Literature overview^{*} | References^{*} |
---|---|---|---|---|---|

Solid particle density | ρ_{s} | kg/m^{3} | 2,650,000 | 2600; 2650; 2340; 2500; 2730; 2650 × 10 ^{9} | [33–38] |

Particles Young's modulus | E_{p} | MPa | 78,000 | 4.61; 26,600; 300;568–2400; 3000–30,000; 4000–75,000; 63.9–120.3; 70,000 | [34, 35, 30, 39–41, 32, 38] |

Particles number | N | – | 46,095 | 2658; 6000; 205–1915; 4150; 10,000; 1 | [33, 42–46] |

Friction coefficient between particles | Μ | – | 0.84 | 0.5; 0.5; 0.25; 0.6; 0.1–0.9; 0.5; 0.7; 0.66; 0.5 | [33, 47, 36, 44, 48, 49, 37, 50, 51] |

Friction coefficient between particles and walls/ porous stone | Μ | – | 0.3 | 0; 0; 1.0; 0.5; 0.5; 0 | [33, 36, 44, 49, 37, 51] |

Rolling friction | μ_{r} | – | 0.04 | 1.6 × 10^{-5}; 0.15; 0.01; 0.05 | [44, 52, 48, 53] |

Time step | Δt | s | 2 × 10^{-7} | 1 × 10^{-7}; 1.25 × 10^{-6}; 2.02 × 10^{-8}; 4.44 × 10 ^{-8}; 1 × 10^{-4} | [33, 44, 54, 55] |

Loading velocity | V | m/s | 4.84 × 10^{-3} | 0.1; 5 × 10^{2}%;0.1 ^{*} × 10^{-4} ε/s | [56, 49] |

All DEM numerical simulation calculations ran on cluster due to faster calculation process. In this case, it is possible to separate modelled oedometer into four quarters where each quarter is calculated by separate computer (**Figure 11**). In DEM simulations, parallel calculations of clusters are widely applied. Using these clusters, it is possible to calculate much bigger simulations in the same period of time, as with a single computer. Nevertheless, applications of clusters in DEM simulations do not provide enough calculation capacity. More detailed explanation of calculations with DEM is presented in Refs. [1, 30].

Used modelled particles shapes in numerical simulations have three shapes: ideal sphere, particle recreated from two spheres and particle recreated from three spheres (**Figure 12**). Recreation of particles shapes is based on experimental testing results obtained from morphology parameters determination part.

Providing DEM numerical simulations with different recreated particles shapes, it is possible to obtain particle shape influence for compression results. The analysis of influence of simulated particles shape on compression results is presented in **Figure 13**.

Analyzing **Figure 13**, it is evident that the compression curve which is obtained from the simulations with ideal sphere particles has a smooth curve. With increasingly complex particle shapes, stress jumps arise in the numerical simulation. These stress jumps appear due to faster particles repositioning than porous stone compression velocity. To increase porous stone velocity is not recommended, because it is possible to have not a static but dynamic modelling effect.

## 5. Comparison of experimental and numerical simulation tests

Reliability and accuracy of numerical DEM modelling of sand soil behaviour depend on the modelled soil and discretization level of particles. In this research work, the same numerical granulometric curve as obtained experimentally was used (**Figure 14**).

The analysis of obtained experimental and numerical morphology parameters is given in **Figure 15**. Comparison of results revealed that to recreate morphology parameters, involvement of different quantity of spheres for single‐particle subscription is necessary: 8–11 spheres for form coefficient, 7–11 spheres for sphericity, 8–12 spheres for circularity, 3–4 spheres for angularity, respectively. Analysing angularity recreation with different quantity of spheres, the angularity coefficient decreases when more than seven spheres are used for single‐particle subscription.

Validation of experimental and numerical (*E _{p}* = 78 GPa and

*ρ*= 2,650,000 kg/m

_{s}^{3}) compression results is presented in

**Figure 16**. One can note that used high modelled particles density increases calculation speed and does not have any effect for contact mechanics due to the first Newton law of motion.

The analysis of compression results revealed (**Figure 16**) that the increased spheres quantity for single‐particle discretization lead the numerical compression curve to the experimental one. For a more qualitative match of experimental and numerical curves, it is better to use the improved experimental compression curve which is closer to the numerical results. In this case, the curves of numerical and improved experimental results get closer to each other. Usually, only numerical tests are corrected with non‐normal coefficients to get the same result curve as in the experiment. In this case, it was improved experimental test to get experimental compression curve more closely to numerical one. This effect was obtained realising the same compression procedure in experiments and in numerical simulations.

## 6. Concluding remarks

This research work presents an experimental investigation of morphology parameters of sand particles, experimental oedometer tests, numerical DEM modelling of oedometer tests, and comparison of experimental and numerical tests. Based on experimental and numerical testing findings, the influence of morphology parameters on mechanical properties of soil is analysed. It is obvious that this research work represents only one soil type (sand). For other soil types as clay, silt, etc., it is not possible to apply obtained results. Nevertheless, the main idea of this research work shows the importance of results validation between experimental and numerical ones. Only comparing results, it is possible to choose the same testing procedure in experimental and numerical simulation. Provided numerical simulations in this research work still are not perfect and needs to be improved. For example, stress jumps, which appear in numerical simulations, are doubtful. The explanation of the stress jumps simple—it is necessary to simulate a bigger particles quantity. In this case, particles for repositioning have not so much space and do not appear big stress jumps. Another way to get smooth and nice numerical simulations compression curve (as in experiment)—simulate particles as ideal spheres. This is what the author wants to show for the readers.

Nowadays simulations oriented on two main things:

Simulate test as simple as just it is possible and accept so far not real particle shape, physical and mechanical properties, but to have smooth compression curve as in experiment. In this case, it would be not possible to use the same input data in DEM program if the modelled soil or soil particles shape would change;

Simulate test as just it is possible with more natural and close to reality parameters. This would make calculations very difficult and calculation time increased. But if the soil or soil particles shape would be changed, it would be possible to run calculations without any change of physical and mechanical properties (just make some changes for modelled particle shape).

In both cases, it is recommended to check the numerical simulations with experimental ones. The influence of the above‐mentioned factors, contributing the accuracy of numerical modelling is a future trend.