The Evaluation of the Parallel Bond ’ s Properties in DEM Modeling of Soils

In this study we conducted the sensitivity analysis of the parallel bond used in the discrete element method (DEM, 3D) based soil model. We researched those parameters which simulate real soil physique attributes. In our investigations we modeled the inhomogeneity by the differentiation in particle size, the soil moisture condition by the parallel bond’s radius R, the cracking tendency by the bonding stiffnesses k k n s , and the air phase by the pore volume of the particle block. We based the validation of the simulation on the results of the simple direct shear box test which were performed in laboratory environment. We analyzed the effect of micromechanical and macromechanical parameters, used in the modeled particle block, with the use of direct shear box and triaxial shear simulations. After the recalculation of micromechanical parameters k k k k n s n s , , , , , , λ σ τ   ( ) we analyzed the effect of the adjustable macromechanical parameters E E c c , ,γ ( ) used in the block with triaxial shear simulation. Based on the comparison of the shear simulations’ results the accuracy of the recalculation of the parallel bond’s micro-macro parameters is proved by the good correlation of the Coulomb failure criterion lines (σ-τ).


Introduction
It is important for the developing engineers of tillage machines to get fast and accurate feedback on the functionality and effectiveness of new constructions, because of the evolved market competition and the increasing needs of users.We can reconstruct field measurements by the investigations conducted in simulated environments.To investigate the micro-and macro cracks appearing in soils we need to use discrete element methods (DEM), instead of using continuous material models (FEM).Nowadays, because of the development of computer science, the DEM method is the most promising approach for modeling the soil-tool (sweep) interaction.
In case of the soil-tool dynamic connection the cracks appearing in the soil affect greatly the size of the emerging deformation zone in front of the tool and thus the tool's towing resistance.
It is hard to set the mechanical parameters of the "synthetic" materials used by the DEM method so as they are almost the same as the mechanical parameters of real materials [2].Whereas the input parameters of continuum material models can be given directly from the results of laboratory measurements (Young's modulus, Poisson ratio, stress, etc.), the macro parameters of the material used in case of the DEM model depend on the features of the micro particles' connections.The input parameters in the DEM model are the parameters of the microscopic particles' connections (for example the spring stiffness between two particles, friction factor, etc.), that is why it is important to calibrate our model with the micro parameters.The parallel bond contact model was used to describe the liquid bridges between the soil particles and to set the so-called capillary effect in the model.The use of the parallel bond contact model enables the realistic cohesion in the DEM model.
The micro parameters of the particle connections , , , , , , λ σ τ   ( ) can't be calculated directly, therefore in the previous DEM models were tried to set the appropriate similarity between the experiments and simulations [3,4,5,6].We can compare the micromechanical parameters of the particle block with the results of the laboratory experiments.The work of [7,8] illustrate well the usage of the material's micromacro parameters validated by modeling.The direct shear simulations used in case of the DEM method can only partially define the macromechanical properties of a given particle block.The triaxial shear simulations were used for stones and rocks, but haven't published results about simulations on soils yet.DEM researches until now haven't focused on the correlations between the soil mechanical parameters measured with direct shear box test and the shear box simulations' micromechanical and triaxial shear simulations' macromechanical setup parameters.
Sadek and his colleagues also used the DEM method in their soil model.Through their investigations they validated their shear test made by the DEM method with laboratory direct shear equipment, paying particular attention to the soil moisture condition [9].
Tsuji and his colleagues modeled the soil moisture condition by simulating the liquid bridges (parallel bond) between the particles and made measurements of the connection of a bulldozer blade and soil [10].
Based on our literature study it can be said that there is no researcher who developed a method, which is verifying the direct shear box and triaxial tests' accuracy of the results to one another, to be able to analyze both quantitatively and qualitatively the sheared particle block.Based on the tuned results of the two tests we can say that we can make the measurable macro and adjustable micro properties in the discrete particle block comparable and convertible.
Through our research we set our goal to iteratively tune in the mechanical parameters of the real soil and of the soil model made by the discrete particle block, throughout this procedure we investigate the correlations between the micro- , ,γ ( ) parameters using the PFC 3D (ITASCATM, USA) software.The goal of the calibration is to compare directly the macromechanical parameters with the real soil mechanical parameters.
Through our work we conducted sensitivity analyses to define the optimal parallel bond model parameters.In the study take account the parameters of parallel bond's radius, bonding stiffness, friction factor between the particle, the porosity and the coord number to investigate the effect that can be consider in the real agricultural soil physics.In addition the use of the mentioned parameters the authors could modelling the soil water content and effect.

Material and method
The shear strength of a given soil is usually defined by laboratory direct shear or triaxial tests.In the former case a so-called shear box is used, this gives the connection between the normal (σ) and shear stress (τ), which describe the Coulomb Law.
In case of triaxial shear the connection can be defined from the Mohr-circles' envelope.From these curves, in both cases the soil's internal friction angle and cohesion can be determined.The triaxial tests are the most reliable, because tests made by using shear boxes are influenced by many parameters, for example by the shape and cross-section of the shear box, the friction between contacting geometries, etc. [11].

Parallel-bond model
Through our research we applied the parallel bond model to analyze the crack propagation and other soil mechanical parameters in soil models (Fig. 1).Capillarity effects can be seen in unsaturated soils, also liquid bridges form between soil particles, which greatly affect cohesion.The parallel bond model connects the discrete particles using a given geometry (circular cross-section, solid beam) [12].Relative movement of connected particles can cause axial and shear forces ) and also torsional moment (M t ) in case of parallel bond connection [13].All in all, parallel bond model can transfer forces and moments also [2]: where F M  The micromechanical parameters of the parallel bond model can be given by the following: the radius multiplier factor λ , which defines the ratio between the parallel bond's cross sectional radius and the radius of the smaller particle in the connection, E c the parallel bond's Young's modulus, k k n s the parallel bond's normal and shear stiffness ratio, the parallel bond's normal and shear directional strength σ c and τ c .We set up the connection between the particles with the help of the contact Young's modulus E c (from which we defined the k n / k s normal and shear stiffness of the particle) and the friction (1) (2) factor (μ) between the particles.In case of the basic connection we had to adjust three micro parameters (E c , k n /k s , μ) and five parameters for the parallel bond λ σ τ , , , , During the validation, the mechanic property of the modeled soil was determined based on the bond strength (micro cohesion) σ τ C C , ( ) used to set the order of magnitude of cohesion within the set.Bond stiffness (Ē c ), used to define the material's tendency for cracking was set relative to the bond strength (Eq.( 4)).It can be proved that the value of Ēc directly affects the dissolution of parallel bonds.Soil compressibility, i.e. the soil's Young modulus (E c ) and Poisson coefficient (γ) can be specified by the value and direction of k n and k s (Eq.( 3)).The magnitude λ ( ) of the bond's radius control the size R ( ) of the liquid bridges between particles the use of Eq. ( 5).
We investigated the convertibility of the macro-and micromechanical parameters of the parallel bond model in order to make the triaxial simulation convertible and comparable to the direct shear box simulation.The setup values can be connected to the geometry of the parallel bond, which is similar to a solid beam with circular cross section going from particle center to particle center.
We can calculate the DEM micromechanical parameters from the macromechanical parameters with a good approximation.It means that from the Young's modulus (E c ) we can define the normal (k n ) stiffness according to Eq. ( 3).
The Young's modulus E c ( ) of the Parallel Bond's define the bond's stiffness k n ( ) in normal direction according to Eq. ( 4).
The parallel bond's normal and shear stiffnesses k n and k s can be expressed by the ratio of the stiffness and the area of the beam's cross section, which creates the connection, as follows [2]: The radius is: The length of the beam is the following [2]:


The shear stiffnesses k s , k s ( ) were set in the proportion of the Poisson ratio (γ).The use of macromechanical parameters can define the micromechanical parameters in the proportion of the the particle radius.

Direct shear test
The validation of the soil mechanical parameters is made possible with the use of the INSTRON 5581 type floor-installed universal mechanical direct shear box test equipment (Fig. 2).We filled the shear box volume with the moist soil used for the testing and then we closed it and applied the preload by a loading mass, which granted the permanent tension load through the measurement.
Because of the displacement of the universal mechanical measurement equipment's cross-head, the top and bottom sides of the measurement box moved relative to each other, while we fixed the force required for the displacement and the value of the displacement itself.From these data the correlation between the normal (σ) and shear stress (τ) can be drawn.

Direct shear box test in DEM
The numerical shear box's plane surfaces are made out of the two half boxes' inner faces, the boundaries associated with them and the top tension plate (Fig. 3).The planes of the shear box are defined by three parameters: the plane normal stiffness, shear stiffness and the friction factor on the walls (Table 1).We created the investigated particle block with gravitational settling method and uniform distribution (Fig. 3 a,b) in the section between the wall elements.We continued creating the elements until the material reached the prescribed porosity.The used shear velocity was the same in case of the measurements and simulations also (300 mm/min).
We investigated the behavior of the particle block using the same particle sizes, thus a homogeneous block (r=3 mm) and using different particle sizes (r=2.55 -3.45 mm) (Table 1).
We analyzed the behavior of the particle block throughout the direct shear box test and direct shear simulations taking into account four different confinements (48700 Pa, 54100 Pa, 66700 Pa and 79400 Pa).We applied the preloads in the direct shear simulations with the servo mechanism that was developed in fish programming language by the authors.The bulk density of the laboratory shear test soil sample was not measured.In the simulation the bulk density was 1560 kg/m 3 [9].
We conducted quality investigations in case of the sheared particle block throughout which we tested the change of the coord number, the value of the sliding fraction and the porosity.The coord number shows us how many particles in general are connecting to a certain one in the particle block.The sliding fraction is defined as the fraction of contacts contained within the measurement sphere that are slipping [2].Slip is assumed to be occurring at a contact if there is no contact bond and the magnitude of the shear component of the contact force is within one-tenth of a percent of the maximum allowable shear force [2].The porosity shows with what volume ratio particles are filling in a defined measuring sphere (given with its center and radius).

Triaxial shear test in DEM
The calibration method of Tamás and his colleagues was used the setting of the DEM model macromechanical parameter [1].We did not conduct a real triaxial test throughout the research.We have taken into account the results of the direct shear box tests in the triaxial tests for which we fit the σ-τ failure envelope using Mohr's circles.The triaxial test made available the convertibility between the micro-macro mechanical parameters.The triaxial shear test is built in function in the PFC 3D .The dimensions of the test specimen throughout the validation: 0.5mx1.0mx0.5m,with 3188 particles (Fig. 4).We defined the minimum coord number (=3) of one element.In the triaxial test the bulk density was 1850 kg/m 3 .
The top and bottom walls were compressed, while the side walls ensured constant preload throughout the triaxial shear simulations thanks to the servo mechanism (Fig. 4).In the triaxial shear simulations in the side walls were taking into account the same confinements were used in the laboratory direct shear box test and in the direct shear simulation (48700 Pa, 54100 Pa, 66700 Pa and 79400 Pa).
In the course of the triaxial test we had to give macromechanical parameters which were calculated from the micromechanical parameters: the contact Young's modulus (E c ), the bonding Young's modulus E c ( ) , the bond's strength σ τ ).
(and the ratio of normal and shear directional components in all cases of the parameters).

The results of investigations
In order to have better knowledge and to be able to use more accurately the parallel bond model to simulate the soil moisture content we performed a sensitivity analysis investigating the effect of the bond's radius, the bond's stiffness and the friction factor between the particles on the quantitative (peak stress) and qualitative (number of cracks) change of the specimen.

The results of the shear simulations
At the initial phase of the direct shear simulation the shear stress (τ) is increasing steeply, almost linearly.This phase is the plastic soil deformation zone.At the end of the first phase the τ stress reaches its maximum (τ B ), this is when the proper shear of the soil cross-section starts.Throughout the progress of the shear process the τ stress gradually decreases until a certain value, then it practically stays constant.After the end of the shear process the soil is sliding on soil, therefore the τ stress is defined by the internal friction angle of the soil.The value of the maximum shear stress:

= + = + τ µσ
Where c -soil cohesion (N/m 2 ), μ -friction factor (soil on soil), τ s -sliding shear stress (N/m 2 ), σ -tension stress (N/m 2 ).The splitting of parallel bonds can be seen on Fig. 5a, which shows well the cracks appearing along the sheared cross section.The velocity vectors of the particles located in the bottom half of the shear box can be seen on Fig. 5b.The top and bottom halves of the shear box had moved 30 mm relative to each other in the tests and simulations also.The results of the triaxial shear test can be seen on Fig. 7a,b.The splitting and disappearance of the parallel bonds from the specimen marked with black can be seen quite well, these outline the crack planes which formed during the shearing.
We conducted a sensitivity analysis on the setup parameters of the parallel bond, which parameters can be traced back to real soil mechanical properties, throughout this we also analyzed qualitative changes in the investigated block.In the course of the sensitivity analysis we investigated the effect of the parallel bond's radius, the bonding stiffness, the friction factor between the particles, the porosity and the average particle bonds.

The effect of the parallel bond's radius on the soil model's mechanical parameters
The moisture capacity in the soil model can be defined by the parallel bond's radius, similarly to the real unsaturated soils.This radius defines the cross section of the liquid bridges connecting the particles.The increase of the radius shows us that in the available range the bonding force between the particles is rising (Fig. 8).We haven't applied values over 1,0 for the radius of the parallel bond, because this kind of arrangement does not occur in reality.In the course of the quantitative investigations we examined the shear peak stresses using different preloads.It can be seen from Figs. 8a,b that the bond force between the particles increases when the radius increases, thus the cohesion in the material rise also in case of homogeneous (Fig. 8a) and inhomogeneous (Fig. 8b) material models.The effect on the coord number of the change of the parallel bond's radius as a function of the preload shows that in case of bigger radii the coord number will be also bigger, this means that a particle will connect to more and more particles (Fig. 9a).The increase of the preload makes this tendency more intensive.The sliding fraction clearly decreases in the smaller preload range when we increase the value of the bond's radius.It seems that in case of bigger preloads, the parallel bond's radius has a lesser effect on the sliding fraction (Fig. 9b).
The porosity clearly decreases due to the increase of the preload; this is not surprising, because we condense the material even more.Because of this the change of the parallel bond's radius shows little effect on the porosity (Fig. 9c).Using smaller preload, we get a saddle surface, thus the porosity has a minimum at the value of 0.5 [-] radius.

The effect of the porosity on the mechanical parameters of the particle block
We analyzed the effect of the pore volume's percent ratio in order to model the air phase in the soil.We investigated the soil model's cohesion and internal friction angle on the range of 30-40% porosity, which showed good correlation with previous research results [15].It can be observed that in both cases of the homogeneous (Fig. 10a) and inhomogeneous (Fig. 10b) material models the increase of the pore volume, that is the decrease of the porosity, the internal cohesion of the soil clearly increases, while it has a lesser effect on the internal friction angle.
In the investigated ranges of porosity we found that the coord number has a minimum value in the porosity of 36-37% (Fig. 11a), at the same time the porosity change has a maximum in this level (Fig. 11a,c).We can observe on the Fig. 11b the change of the sliding fraction.

The effect of the friction factor on the mechanical parameters of the soil model
After the investigation of the effect of the bond's cross section, we have run a sensitivity analysis on the friction factor defined between the particles.By the analysis of the friction factor (µ) between 0.2-1.0values we determined that it has a great effect on the cohesion in the particle block.Based on the Fig. 12 we can declare that the adjusted friction factor between the discrete particles did not have any effect on the internal friction angle.We can see on the Fig. 12 that the change between 0.2-1.0values of the friction factor clearly increases the material's cohesion, while the slope is nearly the same, thus it has a lesser effect on the internal friction angle of the material.
In case of smaller confinements clearly larger values of the coord number belong to smaller values of the friction factor (Fig. 13a).This tendency changes in case of larger confinements, here a saddle surface appears.A local minimum in the coord number arises at around 0.7 friction factor value.
Bigger friction factor clearly decreases the sliding fraction no matter what value the preload is (Fig. 13b).This is in correlation with the increase of the cohesion, because this way fewer bonds disappear between the particles.
A connection is not observable between the friction factor, the preload and the porosity in case of smaller preloads (Fig. 13c).When the preload is large and the friction is small a local minimum value of the porosity clearly appears on the investigated range.

The effect of the parallel bond model's setup parameters on the bond's cracking tendency
The cracks in the soil model mean the disappearance of the parallel bonds (microcracks).Through the quantification of the cracks we were able to investigate the crack sensitivity occurring by the settings made in the DEM model.In case of the triaxial test, we can clearly see the crack, which is built up by microcracks (the disappearance of bonds between particles), running through the specimen.We investigated with the triaxial simulations how the change of the parallel bond's radius affects the block's mechanical parameters.Inside the available range (0.2-1.0) a bigger peak stress (Fig. 14a) and crack number (Fig. 14b) belong to a bigger bonding radius, which gives back the real experimental results.Setting the parallel bond's radius R ( ) to a smaller value models the decrease of the soil moisture condition, in other terms it models the decrease of the liquid bridge's size between two connecting particles.The cohesion only increases to a certain level as a function of the moisture condition, in case of real soils also [16].
Because of the preliminary simulation the parallel bond's stiffness for the crack sensitivity analysis in the range of Ē c = 0.8e6-6e6 Pa/m was set in accordance with the parallel bond's strength σ τ ).In the course of the simulation, for a given stress strength if we set the bonding stiffness higher than the materials cracking tendency gets higher (Fig. 15a) and we obtain higher peak stresses (Fig. 15b).We can define the development of the cracks in the specimen from a gel-like state (no cracks) to a crumbling state (all connections disappear).
The number of cracks, because of the effect of the preload's and the friction factor's change, was the largest when we investigated the particle block with the largest preload and smallest friction factor.This can be explained by the cohesion increasing effect of the friction factor.
A similar tendency can be observed for the peak stress (Fig. 16b).In this case the value of the peak stress steeply rises as the effect of the increasing preload and friction factor, until it reaches a limit, then the fitted surface almost turns into a plane.Based on this it can be said that over a certain preload and friction factor the peak stress does not increase.

The effect of the coord number on the parallel bond's cracking tendency
In the simulated triaxial shear tests, we investigated the changes in the particle block using the so-called number of the particles' average connections (coord number), which has a similar effect as the porosity.We investigated the effect of the average particle connections on the results of the triaxial simulations in the range of 1.00-5.00particle connections.
A peak (coord number=4) can be observed after a mild rise when we increase the coord number in the range of 1-3, then a decreasing tendency can be observed (Fig. 17a).The peak stress increases with the increase of the particle connections (Fig. 17b), which result reflects well previous research results [17,18].

The comparison of the results of the real and simulated shear tests
We can calculate the DEM model micromechanical properties from the macromechanical parameter with a good approximation.In this case the the shearing simulations are comparatively (Fig. 18).The comparison of the results of the triaxial shear simulation, the direct shear simulation and the direct shear test can be seen on the Fig. 18.The good correlation of the failure lines (σ-τ) verifies the convertibility of the microand macromechanical initial setup parameters.In the course of the investigation, we used the triaxial shear simulation's macro parameters in the direct simulation also.Nevertheless, we have to mention that the internal friction angle showed a higher value of 1-3 degrees.
In case of the validation of the triaxial tests we started from the bond's strength σ τ C C , ( ) between the particles (microcohe- sion) to setup the mechanical parameters of the soil model, with which we were able to adjust the order of magnitude of the cohesion.We set the bonding stiffness ( E c = 2e5 Pa) according to the stress strength between the particles, we defined the crack tendency this way.It is provable that this value has a direct effect (pb_kn=5e6 Pa/m, pb_ks=2,5e6 Pa/m calculated from E c ) on the disappearance of the parallel bonds, thus on the development of crack planes.The compressibility of the soil, the value of the Young's modulus (E c =1e6 Pa) and the Poisson ratio can be given by the ratio and value of the k n and k s stiffnesses (stiffnesses of springs between two particles).The adjusted value of the stiffnesses k n =3,14e4 N/m, k s =1,57e4 N/m at the contact of the particles has returned inside the given order of magnitude the behavior of sandy soil from the soil trough using particles with 20-26 mm radii.The beforehand mentioned particle radius dependent micromechanical parameters expressed by the macromechanical parameters independent from the particle radius: Young' modulus (E c ): 1e6 Pa, Poisson ratio (γ): 0.281.

Conclusion
In our research, we investigated the availability of discrete element modeling of agricultural soils used in the soil-sweep interaction model.Based on the conducted simulations it can be said that the DEM method with parallel bond settings is able to model real soils with high accuracy.We modeled the unsaturated agricultural soil's moisture condition with the radius of the parallel bond and the air phase with the pore volume.We studied with further sensitivity analyses the effect of the friction factor and the bonding stiffness on the cohesion and crack tendency of the particle block.
We based the validation of the simulations on the laboratory direct shear box simulations.We analyzed the effect of the adjusted micromechanical parameters in the particle block with the virtual direct shear box simulation.After the recalculation of the micromechanical parameters k n ,k , , , , , s R k k we analyzed the effect of the adjusted macromechanical parameters in the block using the virtual triaxial simulation.The accuracy of the recalculation is verified by the good correlation of the failure lines (σ-τ) which were based on the results of the shear simulations.The result of our research grants the application of the macromechanical parameters in case of using the parallel bond model in discrete element simulations.We solved the interoperability between the micro and macro parameters.The convertibility grants the utilization of the previous research results of related literatures in the DEM model.

Fig. 3
Fig. 3 a) Block generation by settling, b) the dimensions of the direct shear box used in the validation (200mmx200mmx(2x20mm)

Fig. 4
Fig. 4 Dimensions of the test specimen used for the triaxial simulations: 500 mm x 1000 mm x 500 mm, 3188 particles

Fig. 5 a
Fig. 5 a) Splitting of parallel bonds b) velocity vectors while pulling the bottom half box

( 7 )Fig. 7 a
Fig. 7 a) Parallel bonds (marked with black) splitted during shearing and the shear curve in case of 54156 Pa confinement stress and b) the velocity vectors based on the results of the virtual triaxial shear test

Fig. 8
Fig. 8 Results of direct shearing using different parallel bond cross sections with a) homogeneous and b) inhomogeneous particle sizes

Fig. 9
Fig. 9 The effect of the parallel bond's radius on a) the coord number, b) the sliding fraction and c) the porosity in case of different confinements applied to the inhomogeneous material direct shear simulations

Fig. 10
Fig. 10 The effect of porosity in the particle block on the results of the direct shear box simulations in case of a) homogeneous and b) inhomogeneous material models

Fig. 11
Fig. 11 The effect of the porosity in the investigated block on a) the coord number, b) the sliding fraction, c) the porosity in case of the applied confinements in the direct shear simulations

Fig. 12
Fig.12 The effect of the friction factor (0,2-1,0) in case of the inhomogeneous particle block

Fig. 13
Fig. 13 The effect of the adjusted friction factor in the investigated block on a) the coord number, b) the sliding fraction, c) the porosity in case of the applied confinements in the direct shear simulations

Fig. 14 Fig. 15
Fig.14 The effect of a) the size of the parallel bond's radius and b) The bonding stiffness on the crack sensitivity used in the connections of the particles in case of 48698Pa, 54156Pa, 66742Pa, 79352Pa confinements

Fig. 17
Fig.17 The effect of the coord number (1,00-5,00) on a) the number of cracks and b) the peak stress

Table 1
The parameters of the shear box and triaxial shear test simulations