Analytical and Optimization Study of Propagation and Reduction of Wave in Granular Sand by DEM Method

Amin Moslemi Petrudi^* and Masoud Rahmani

Department of Mechanical Engineering, IHU University, Tehran, Iran

E-mail: amin.moslemi2020@gmail.com; msrahmani@ihu.ac.ir

^*Corresponding Author

Received 20 January 2021; Accepted 10 March 2021; Publication 06 May 2021

Abstract

In this research, the discrete element method has been used to analyze wave propagation and to investigate the factors affecting wave reduction in granular soils. The method of discrete elements is important because of the possibility of preparing completely similar specimens and examining the effect of changes in a certain parameter on the Behavior of the specimens. This method also provides an understanding of the changes that have occurred at the micro-scale of granular materials that are not achievable with other laboratory and numerical methods. To model the specimens, a set of disks with specific granulation has been used for two-dimensional studies. PFC 2D software has been used to perform simulations and related analyzes such as interparticle force. The DEM code in MATLAB is used to check the wave depreciation. In this research, the optimization process was performed using experimental data and the Taguchi method using the DEM method. The results of this study show that there is a direct relationship between the number of particle set contacts and the wave propagation speed. Also, material properties such as particle density are the most important parameters affecting wave velocity. The results of the method (DEM) are done with PFC 2D software and a comparison between the results of this method with the solution methods used by other researchers is shown to be a good match.

Keywords: Discrete element method, wave, granular, propagation.

1 Introduction

Granular materials are composed of distinct particles that exhibit complex macroscopic behavior against external loads. Soils are also materials composed of particles of different sizes and their behavior is determined by the forces between these particles. However, this feature is not usually considered in the modeling. Soil interparticle forces include forces due to boundary conditions, interparticle forces (contact forces) that the relative balance between these forces reveals different aspects of soil behavior. The phenomenon of wave propagation plays an essential role in various dynamic issues such as soil-structural seismic interaction, liquefaction, and vibration. Understanding the local effects of the site on strong ground movements and assessing the response and deformation of the ground to strong movements is critical for critical structures and facilities. The study of wave propagation in granular materials also has important industrial applications. Granular materials are used to absorb shock waves during the transmission of heavy equipment and to isolate sensitive equipment from earthquakes. They are also used in the manufacture of ceramic components that require dynamic compaction of ceramic powders. In all these applications, it is necessary to study the wave velocity and the nature of its propagation in granular materials. Research on pressure wave transmission in granular soils has been performed by various researchers. A. Levy [1] presented his studies in an article entitled Shock Wave Contact with Granular Materials. In this paper, the state model developed for porous materials was used for granular materials, and the one-dimensional state of these equations by The numerical code is solved and shows that it is in good agreement with the experimental work of other researchers. A-Britan et al. [2] studied the reduction of the shock wave in grain filters. Because most grain filters made of gravel and sand are considered spherical, in this study the effect of parameters such as particle shape has been neglected The presence of an air gap between the granule filter and the protected surface prevents direct contact between the particles forming the grain filter and the protected surface. Increasing the filter length or decreasing the particle diameter reduces the amplitude of the output shock wave, the particle density, and lateral friction do not have much effect on reducing the gas pressure in the filter. Vitali F. Nesterenko [3] in a laboratory work entitled Shock Wave Depreciation Using Soft Compressed Materials investigated the effect of using excavation in a metal tank and observed that using this material greatly reduced the deformation of the tank and its degradation. It is prevented. Another researcher in this field, Stephen R. Hostler [4], in his dissertation has investigated how wave propagation in granular materials using numerical and experimental methods. He performed the simulation in two dimensions and used the Discrete Element Method (DEM). In this study, the effect of substrate vibration on wave propagation was investigated. One of the results of their work was that with the vibration of the substrate, the velocity of the passing wave generally increases. Susana Serna et al. [5] also studied the reduction of the blast wave in the grain flow, they studied the hydrodynamics of the blast wave in the grain gases. C.M. Wensrich et al. [6] investigated the shock wave in granular material. He also used DEM simulation in his work. The mechanism of shock formation was observed in a series of DEM simulations using the Hertz-Mindlin contract law. This model with a linear contact model between particles concluded that the use of the linear contact model for dynamic loads such as impact and shock in granular materials is not appropriate. C. Guéders et al. [7] have simulated the explosion wave depletion in a bed of alumina ceramic beads by considering pressure and impulse. Numerical simulations have been performed by Zamani et al. [8] to study shear wave propagation in dry grain soils. In this study, periodic boundaries are used to simulate an infinite system on both sides. Such boundaries allow particles to reach an imaginary adjacent surface from the parallel region and simulate an infinite periodic system in lateral directions. In this study, they used the DEM method for simulation. Ngoc-Son Nguyen et al. [9] Numerically studied the shock wave dynamics in a chain of granular materials. Comparison of their work with experimental results showed that the multiple impact model is well able to describe the shock wave dynamics in the grain chain. W.D. Neal et al. [10] in a study investigated the possibility of sustained shock waves of brittle particles such as sand. They examined the effect of particle dispersion on the shock wave velocity by varying the grain bed thickness and input stress. They also showed that the wave velocity increased with the increase of the applied stress and also showed that the wave passage time is proportional to the thickness of the granular substrate. Nesterenko [11] in his book entitled Dynamics of Heterogeneous Materials deals with topics such as wave propagation in the granular matter in two and three-dimensional states, wave propagation in a chain of grains which is a simplified form of granular material to study wave dynamics. He has also propagated shock waves in porous and multilayer materials. In this book, he has used both experimental and numerical methods. A.-N. Rotariu et al. [12] investigated the measurement of pulse attenuation in environments of porous granular materials by the experimental test method. Permanent Sheets Estimate Wave Depreciation. K.Kandan et al. [13] investigated the instability of the grain bed surface in shock loading. In their work, the onset and growth of instability in granular materials loaded by air shock waves were investigated through shock tube tests and numerical calculations. In their study, the Lagrangian-Eulerian solver and Drucker-Prager model for granular materials have been used. Liang et al. [14] studied the shock waves in a one-dimensional bead produced by a constant velocity impact in a short period using numerical simulation. In their study, a cylinder and piston were used to generate the wave. Rahmani and Moslemi [15] investigated the shock wave attenuation in sawdust and mineral aggregate aggregates by simulated discrete particle hydrodynamics (SPH) and performed an experimental test. In their study, a sandwich panel with an aggregate core using a shock tube Exposed to the blast wave.

In this research, the discrete element method has been used to analyze wave propagation and to investigate the factors affecting wave velocity in granular soils. The method of discrete elements is important because of the possibility of preparing completely similar samples and examining the effect of changes in a certain parameter on the behavior of the samples. It also provides an understanding of micro-scale changes in aggregates that cannot be achieved with other laboratory and numerical methods. For modeling purposes, samples were created from a set of discs with specific grain sizes for two-dimensional studies. PFC 2D software has been used to perform simulations and related analyzes. Finally, the simulation performed by DEM should be compared with the results of other laboratory experiments to measure the accuracy of numerical modeling. To do this, in the present study, the results of pressure tests [15] have been used.

2 Discrete Elements Method

The discrete element method is one of the numerical methods that has been developed specifically to model the behavior of discontinuous systems such as grain environments. This method was first used by Cundall (1971) to model the progressive failure of rock slopes. Cundall & Strack (1979a) then developed their original method for modeling a two-dimensional set of circular plates [16]. Ghaboussi & Barbosa (1990) developed a method for modeling the behavior of granular materials and various geometric shapes [17]. In this method, each particle is modeled physically from a set of particles separately, taking into account the geometric conditions of its surface and a description of its physical state (location, direction, volumetric forces, etc.). The motion of particles is determined by the forces applied to them and by Newton’s second law. The physical and mechanical conditions of particles are determined to mathematically model them using particle micromechanical relationships.

3 Equations of Motion

The law of transient motion for a single particle is written as follows:

{\ddot{u}}_{i} + α {\dot{u}}_{i} = \frac{f_{i}}{m} + g_{i}

(1)

In the above relation ${\ddot{u}}_{i}$ is the acceleration of the motion of the center of the particle, $α$ is the viscous energy absorption constant, $f_{i}$ is the sum of the forces acting on the particle due to the contacts, m is the mass of the particle and $g_{i}$ is the acceleration vector of the earth’s gravity. Similarly, for the rotational motion:

{\ddot{θ}}_{i} + α {\dot{θ}}_{i} = \frac{M_{i}}{I}

(2)

And in this relation ${\ddot{θ}}_{i}$ is the angular acceleration of the center of the particle and $M_{i}$ is the total anchor on the particle.

Using the central finite difference method and integrating the equations of motion, the relative velocity and angular values in the middle of each time interval are obtained:

${\dot{u}}_{i}^{t}$	$= \frac{1}{2} {{\dot{u}}_{i}^{(t - \frac{Δ t}{2})} + {\dot{u}}_{i}^{(t + \frac{Δ t}{2})}}$	(3)
${\dot{θ}}_{i}^{t}$	$= \frac{1}{2} {{\dot{θ}}_{i}^{(t - \frac{Δ t}{2})} + {\dot{θ}}_{i}^{(t + \frac{Δ t}{2})}}$	(4)

For accelerate:

${\ddot{u}}_{i}^{t}$	$= \frac{1}{Δ t} {{\ddot{u}}_{i}^{(t - \frac{Δ t}{2})} - {\ddot{u}}_{i}^{(t + \frac{Δ t}{2})}}$	(5)
${\ddot{θ}}_{i}^{t}$	$= \frac{1}{Δ t} {{\ddot{θ}}_{i}^{(t - \frac{Δ t}{2})} - {\ddot{θ}}_{i}^{(t + \frac{Δ t}{2})}}$	(6)

By placing the above values in the equations of motion (displacement) and motion (angular) and solving the equations for the velocity at time $t + \frac{Δ t}{2}$

${\dot{u}}_{i}^{(t + \frac{Δ t}{2})}$	$= {D_{1} {\dot{u}}_{i}^{(t - \frac{Δ t}{2})} + (\frac{f_{i}^{(t)}}{m} + g_{i}) Δ t} D_{2}$	(7)
${\dot{θ}}_{i}^{(t + \frac{Δ t}{2})}$	$= {D_{1} {\dot{θ}}_{i}^{(t - \frac{Δ t}{2})} + (\frac{M_{i}^{(t)}}{I}) Δ t} D_{2}$	(8)

In the above relations, the coefficients $D_{2}$ and $D_{2}$ :

$D_{1}$	$= 1 - α \frac{Δ t}{2}$	(9)
$D_{2}$	$= \frac{1}{1 + α \frac{Δ t}{2}}$	(10)

The amount of relative and angular displacement changes is calculated according to the following equations:

$Δ u_{i}$	$= {\dot{u}}_{i}^{(t + \frac{Δ t}{2})} Δ t$	(11)
$Δ θ_{i}$	$= {\dot{θ}}_{i}^{(t + \frac{Δ t}{2})} Δ t$	(12)

Finally, the location of the $i$ particle changes to the following value:

u_{i}^{(t + Δ t)} = u_{i}^{(t)} + Δ u_{i}

(13)

It should be noted that the sum of forces and anchors on each particle ( $f_{i}$ and $M_{i}$ ) changes to zero after updating the particle location in each cycle. Then, the above set of steps is repeated for the new location of each particle in the force-displacement relationship.

4 Validation and Modeling

The first step in modeling is to create a specimen with the desired porosity and dimensions. There are several methods for producing a set of particles. The method used in this modeling is the particle production method with increasing radius. In this method, first, in a specific space, several circular particles that will not come into contact with each other are randomly generated. The particles will then be large enough to achieve the desired porosity. Walls are used to creating boundary conditions. Initially, the walls are created as borders so that a specimen with specific dimensions can be created. The initial conditions of stresses in a set of particles depend on how the particles are produced and compacted. To achieve the desired initial stresses, a small amount of change in porosity can be made, because small changes in porosity will cause large changes in stresses. The behavior of the materials is simulated using the contact model in each contact. Depending on the contact model, there may be parallel adhesion and damping. As shown in the figure below, there are three components to defining force-displacement behavior. Figure 1 shows the calculation cycle of the DEM method.

Figure 1 Calculation cycle in DEM method [18].

Behavioral components used in the models include stiffness, slip, and adhesion. Contact stiffness Relate contact forces to relative displacements using the following equation. Normal stiffness, total normal force, is related to total displacement.

F^{n} = K^{n} U^{n}

(14)

Tangential shear stiffness, shear force development are related to shear displacement development.

F^{s} = - K^{s} U^{s}

(15)

Contact stiffness is different for different contact models.

A slip occurs in a relationship between shear and normal force on both sides of the contact, which causes the two particles to slide relative to each other.

Slip behavior is defined by the coefficient of friction ( $μ$ ) at the point of contact. Where $μ$ is the smallest coefficient of friction of two particles in contact with each other.

F_{\max}^{s} = μ | F_{i}^{n} |

(16)

If $| F_{i}^{s} | > F_{m a x}^{s}$ , the slip will occur.

In this study, according to the type of soil studied (granular), adhesion was not used.

The Hertz model of slip and stiffness (which changes as a function of the elastic constants of the two sides of the contact) provides normal forces and overlap.

The Hertz-Mindlenn contact model is a nonlinear contact model based on the theory of Mindlin and Deresiewics (1953) proposed by Cundall (1988). This model is only applicable to circular particles.

The model is defined by two parameters, shear modulus and Poisson’s ratio of two particles in contact with each other. For the Hertz contact model, two parameters Kn and Ks are defined.

The difficulty of normal sequential contact is calculated by the following equation:

K^{n} = (\frac{2 G \sqrt{2 R}}{3 (1 - υ)}) \sqrt{U^{n}}

(17)

The tangential shear contact stiffness is obtained from the following equation:

k^{s} = (\frac{2 {({(G)}^{2} 3 (1 - υ) R)}^{1 / 3}}{2 - υ}) {| F_{i}^{n} |}^{1 / 3}

(18)

Where $U^{n}$ is the overlap value and $| F_{i}^{n} |$ the magnitude of the normal contact force The coefficients of the above two equations are a function of the geometry and material properties of the particles on both sides of the contact. For butter-to-butter contact, the coefficients are calculated as follows:

$R$	$= \frac{2 R^{A} \times R^{B}}{R^{A} + R^{B}}$	(19)
$υ$	$= \frac{1}{2} (υ^{A} + υ^{B})$	(20)
$G$	$= \frac{1}{2} (G^{A} + G^{B})$	(21)

For circle contact with the boundary, the coefficients will be as follows:

υ = υ^{b a l l}, G = G^{b a l l}, R = R^{b a l l}

(22)

In the Hertz model, the sequential normal stiffness is related to the tangential normal stiffness by the following formula [8].

k^{n} \equiv \frac{d F^{n}}{d U^{n}} = \frac{3}{2} K^{n}

(23)

Choosing the right values of material parameters is difficult. For continua codes such as finite elements and finite differences, the input properties of materials such as modules and resistors can be derived directly from the results of laboratory measurements. But this method requires the properties of micro-materials. In this study, a large number of micro-properties that were available were assigned to materials. The rest of the properties are done by modeling and the actual results and achieving the desired results are specified and assigned to the materials. Table 1. Provides simulation parameters.

Table 1 Simulation parameters

Specimen height	40 cm
Specimen width	25 cm
Particle mass $(\frac{kg}{m^{3}})$	2500
Poisson’s ratio	0.22
Modulus of elasticity (GPa)	70.3
friction coefficient ( $μ$ )	0.1
Number of particles	4000

The energy applied to the particles is dissipated through frictional slip. However, this rate of energy consumption may not be sufficient to reach a steady-state in several logic analysis cycles. Therefore, it is necessary to use mechanical damping (local damping, hysterical damping, and viscous damping) to consume kinetic energy. Local damping acts on each particle, while viscous damping applies to each contact. In local damping, a damping force of magnitude proportional to the unbalancing force is applied to each particle. In viscous damping, normal and shear dampers are added per contact. These dampers operate in parallel with the existing contact model and provide forces that are proportional to the relative velocity difference between the two sides of the contact (particle-particle or particle-wall). The damping force parameter is added to the motion equation as follows.

$F_{(i)} + F_{(i)}^{d}$	$= M_{(i)} A_{(i)}; i = 1 \dots 6$	(24)
$M_{(i)} A_{(i)}$	$= {\begin{matrix} m \ddot{x}, & for i = 1 \dots 3 \\ I {\dot{ω}}_{(i - 3)}, & for i = 1 \dots 6 \end{matrix}}$	(25)

Where $F_{i}$ , $M_{i}$ , and $A_{i}$ are the components of force, mass, and acceleration, respectively, and $F_{i}^{d}$ is the damping force.

$F_{(i)}^{d}$	$= - α \| F_{(i)} \| s i g n (v_{(i)}); i = 1 = \dots 6$	(26)
$s i g n (y)$	$= {\begin{matrix} + 1, & if \dots y > 0 \\ - 1, & if \dots y < 0 \\ 0, & if \dots y = 0 \end{matrix}}$	(27)
$v_{(i)}$	$= {\begin{matrix} {\dot{x}}_{(i)}, & for i = 1 \dots 3 \\ w_{(i - 3)} & for i = 4 \dots 6 \end{matrix}}$	(28)

The damping force is controlled by the damping constant ( $α$ ), which can be assigned to each particle separately. This form of attenuation has the following advantages:

1. Only the acceleration parameter is damped. Therefore, there will be no false damping forces due to the steady-state.

2. The damping constant ( $α$ ) is dimensionless.

3. Since the damping is independent of the frequency, different parts of the set are damped equally with different natural periods with a damping constant.

For dense sets, local damping is the best option to achieve equilibrium and perform quasi-static simulations.

The loading method in this simulation is using the movement of borders. The velocity of the wall is based on a sinusoidal function through which a pressure wave is transmitted to the particles attached to the wall, and in the same way, these particles transmit the stress wave to the particles attached to them. Modeling geometry and input wave are seen in Figures and. The input wave vibration frequency in the simulations is 20 Hz and the input wave amplitude is equal to 0.03 m/s. Figure 2 shows how to apply a load to a set of particles and Figure 3 shows the input wave in the simulation.

Figure 2 Apply a load to a set of particles.

Figure 3 Input wave in simulation.

4.1 Validation

The DEM method can also be examined in the form of analytical coding. To perform the simulation validation, we have used the DEM code using MATLAB software. In the DEM model, the same sand properties are considered, stress diagrams. The strain due to axial pressure is obtained by the DEM method and compared with the diagrams obtained from the simulation. Figure 4 shows the stages of particle compaction. The particles are randomly generated in a range of 1 by 2 cm, then they are compressed. A compression diagram is extracted and from it, according to the fixed cross-sectional area, the strain stress is obtained. Particle diameters between 1 and 0.15 mm are random and the number of particles in this section is 75.

Figure 4 Particle compression steps by DEM method.

First, the force-displacement diagram (stress-strain) obtained from axial loading with Dem code is compared with the experimental test [15] shown in the figure. The graphs are obtained from a press with a load of 1.5 tons with a compression speed of 2 mm per minute. Figure 5 shows a comparison of the DEM code with the experimental test [15] and Table 2: shows the specifications of the test material.

Figure 5 Calculation cycle in DEM method [18].

Figure diagram with the experimental result obtained, you can see that the behavior of the diagram is similar to the diagram obtained from the code of the DEM method.

Table 2 Specifications of the test material

Particle Shape	Size (mm)	$C_{0}$ (m/s)	$υ$	$E (G P a)$	$ρ (\frac{k g}{m^{3}})$
Spherical	2–8	5300	0.22	70.3	2.47

5 Result and Discussion

First, the wave transmission in granular soils will be examined, then the factors that may affect the wave velocity will be investigated. These factors include soil porosity, particle surface hardness, density, PDI, and soil granulation, the effect of each of the above factors will be investigated by forming different sets of particles with different parameters. Figure 6 shows the contact forces at times 1, 2, 3, and 4 milliseconds, and Figure 7 shows the velocity wave at different distances from the center at the lower level of the specimen.

Figure 6 Showing contact forces at times of 1, 2, 3, and 4 milliseconds.

Figure 7 Velocity wave at different distances from the center at the bottom of the specimen.

The results of the obtained velocities for the specimens with different porosities are shown in Figure 9. As can be seen, the wave propagation speed increases with decreasing soil porosity. In specimens with less porosity, the number of contacts between the particles increases, which in turn causes the specimen to become stiffer. As a result, this increase in specimen hardness will increase the wave velocity. Figure 8 shows the velocity changes with porosity.

Figure 8 Speed changes with porosity.

This result is consistent with the results of experiments and simulations performed by Zamani et al. That three soil specimens with different porosities were tested under shear waves.

11 Graphs of porosity changes versus time of wave application for different soils are shown. Examining the above diagrams, it is found that by applying a wave to the particle set, the specimen begins to compact, which is faster in high-porosity soils than in low-porosity soils. The reason for this can be due to the presence of more voids in soils with high porosity. Also, due to the smaller number of contacts in the specimen, the particles have more freedom to move, thus having more potential to be compressed. Figure 9 shows the change in porosity with time.

Figure 9 Change of porosity coefficient with time.

The coefficient of friction is another factor affecting the speed of wave propagation. In this section, three soil specimens with a friction coefficient of 0.1, 0.3, 0.5 are made and compared.

Examining the above diagrams, it is found that by applying a wave to the particle set, the specimen begins to compact, which is slower in compacted soils with a higher coefficient of friction. The reason for this can be attributed to the greater freedom of particles with less friction to move on each other and slip. As a result, a set with a lower particle friction coefficient will have more potential for compression.

The velocity-time diagram of the particle is shown in each case of the coefficients of friction. According to the above figures, the obtained wave velocities are equal to 0.158, 0.163, 0.169, and 0.174 for friction coefficients of 0.1, 0.3, 0.5, and 0.7, respectively. Be. According to the presented results, the wave propagation speed increases with increasing the coefficient of friction (particle surface hardness). Figure 10 shows the speed changes with the coefficient of friction.

Figure 10 Shows speed changes with the coefficient of friction.

Structural characteristics of particles are one of the factors affecting the speed of wave propagation. One of the most important of these characteristics is density. In this section, three types of densities of 2400, 2000, and 1800 (kg/m $^{3}$ ) have been studied.

The results related to the wave propagation velocity obtained in specimens with different densities are shown in Figure 12. As can be seen, the velocity of the wave propagation increases with decreasing density. Increasing the density increases the inertia or inertia of the grains, so the particles are less inclined to change their position. Eventually, this increase in particle density will reduce the wave velocity in the particle set. Figure 11 shows the velocity changes with density.

Figure 11 Shows velocity changes with density.

Figure 12 Shows speed changes with PDI.

5.1 Definition of Grain Non-uniformity Coefficient (PDI)

PDI means measuring the amount of non-uniformity of grain or particle size in a specimen. If all the particles have the same size, shape, and mass (Figure 6), the specimen is called monodisperse. If the set contains particles of different shapes, masses, and sizes (Figure 6 the specimen is called non-uniform or polydisperse.

The PDI (polydispersity index) is determined using Equation (29).

$P D I$	$= \frac{M_{w}}{M_{n}}$	(29)
$M_{n}$	$= \frac{\sum_{i = 1}^{n} N_{i} M_{i}}{\sum_{i = 1}^{n} N_{i}}$	(30)
$M_{w}$	$= \frac{\sum_{i = 1}^{n} N_{i} M_{i}^{2}}{\sum_{i = 1}^{n} N_{i} M_{i}}$	(31)

Figure 13 Three-dimensional graph of the ratio of DEM to the experimental process.

Figure 14 Interaction ratio of DEM to the experimental process.

To investigate the effect of PDI on wave velocity, four soil types with 5%, 10%, 20%, and 30% PDIs were studied. The number of particles in the made specimens is equal to 8000.

According to studies, the average number of calls in all four cases is approximately equal to 3.55. Also, the contact force between the particles is almost the same in all four cases.

Table 3 Points predicted by Taguchi method in optimization

Predicted
Response	Mean	Median $^{1}$	Observed	Std Dev	SE Mean

CI for	Mean		99% of		Population

95% CI Low	95% CI High		95% TI Low		95% TI High

Factor	Name	Level	Low Level		High Level

A	A	-1.35484	-1.35484		-2284.65
B	B	-1.75	-1.75		-2291.29
C	C	-2.20161	-2.20161		-1988.58
D	D	-2.48387	-2.48387		-2713.17
E	E	-2.70968	-2.70968		-3331.05
F	F	-2.99194	-2.99194		-3845.53
G	G	-3.27419	-3.27419		-4773.59
H	H	-3.44355	-3.44355		-5289.73
J	J	-3.66935	-3.66935		-5810.03
K	K	-3.83871	-3.83871		-6126.03
L	L	-4.23387	-4.23387		-6438.7
M	M	-4.40323	-4.40323		-7054.92
N	N	-4.74194	-4.74194		-7778.68
O	O	-4.96774	-4.96774		-7682.76
P	P	-5.30645	-5.30645		-7584.35
Q	Q	-5.70161	-5.70161		-8206.37
R	R	-5.98387	-5.98387		-8722.52
S	S	-6.32258	-6.32258		-9136.09
T	T	-6.77419	-6.77419		-9549.67
U	U	-7.05645	-7.05645		-10273.4
V	V	-7.33871	-7.33871		-10585.3

Table 4 Factors considered by the Taguchi method

Study Type	Factorial		Subtype	Randomized

Design Type	Taguchi OA		Runs	64

Design Model	Main Effects		Blocks	No Blocks	Build Time (ms)

Factor	Name	Type	Subtype	Minimum	Maximum

A	A	Categorical	Nominal	–1.35484	–2284.65
B	B	Categorical	Nominal	–1.75	–2291.29
C	C	Categorical	Nominal	–2.20161	–1988.58
D	D	Categorical	Nominal	–2.48387	–2713.17
E	E	Categorical	Nominal	–2.70968	–3331.05
F	F	Categorical	Nominal	–2.99194	–3845.53
G	G	Categorical	Nominal	–3.27419	–4773.59
H	H	Categorical	Nominal	–3.44355	–5289.73
J	J	Categorical	Nominal	–3.66935	–5810.03
K	K	Categorical	Nominal	–3.83871	–6126.03
L	L	Categorical	Nominal	–4.23387	–6438.7
M	M	Categorical	Nominal	–4.40323	–7054.92
N	N	Categorical	Nominal	–4.74194	–7778.68
O	O	Categorical	Nominal	–4.96774	–7682.76
P	P	Categorical	Nominal	–5.30645	–7584.35
Q	Q	Categorical	Nominal	–5.70161	–8206.37
R	R	Categorical	Nominal	–5.98387	–8722.52
S	S	Categorical	Nominal	–6.32258	–9136.09
T	T	Categorical	Nominal	–6.77419	–9549.67
U	U	Categorical	Nominal	–7.05645	–10273.4
V	V	Categorical	Nominal	–7.33871	–10585.3

The presence of particles with different diameters in a set alone cannot change the wave velocity. Changes in particle size can cause changes in wave velocity if they increase the stability of the set (increase the number of medium contacts and make the set of particles more difficult). Figure 12 shows the speed changes with PDI.

The results of evaluating the effect of this parameter are shown in the figure. As can be seen, no clear correlation can be found between PDI and wave velocity. This may be due to other factors affecting the wave velocity.

5.2 Optimization

To achieve the desired results in 21 stages and with different values, 21 data factors are presented. To compare the process of optimizing the experimental results with the DEM method, its data were extracted using the Taguchi method. Figure 6 shows a three-dimensional graph of the ratio of DEM to the experimental process and Figure 7 shows the interaction ratio of DEM to the experimental process. The purpose of optimizing the experimental method and the DEM method is to obtain the optimal points and their best performance. This optimization has obtained appropriate coefficients by determining the objective functions and experimental data. Table 3 shows the points predicted by the Taguchi method in optimization and Table 4 shows the factors considered by the Taguchi method.

6 Conclusion

The results of the obtained velocities were shown for the samples with different porosities. As can be seen, the wave propagation rate increases with decreasing soil porosity. According to the results of wave propagation analysis in soils with different porosity, in samples with less porosity, the number of contacts between particles increases, which ultimately causes the sample to become stiffer. As a result, this increase in sample hardness will increase the wave velocity. According to the study of the effect of the coefficient of friction on wave velocity, it was observed that with increasing the coefficient of friction, the wave propagation velocity decreases. As the density decreases, the wave propagation velocity increases. Increasing the density increases the inertia or inertia of the grains, so the particles are less inclined to change their position. You have also seen that both DEM methods have good ability in modeling granular materials. The results of the method (DEM) are done with PFC 2D software and a comparison between the results of this method with the solution methods used by other researchers is shown to be a good match.

References

[1] Levy, A. “Shock waves interaction with granular materials.” Powder technology 103, no. 3 (1999): 212–219.

[2] Britan, A., G. Ben-Dor, O. Igra, and H. Shapiro. “Shock waves attenuation by granular filters.” International Journal of Multiphase Flow 27, no. 4 (2001): 617–634.

[3] Nesterenko, Vitali F. “Shock (blast) mitigation by “soft” condensed matter.” arXiv preprint cond-mat/0303332 (2003).

[4] Hostler, Stephen R., and Christopher E. Brennen. “Pressure wave propagation in a granular bed.” Physical review E 72, no. 3 (2005): 031303.

[5] Serna, Susana, and Antonio Marquina. “Capturing blast waves in granular flow.” Computers & Fluids 36, no. 8 (2007): 1364–1372.

[6] Wensrich, C. M., and R. E. Stratton. “Shock waves in granular materials: Discrete and continuum comparisons.” Powder Technology 210, No. 3 (2011): 288–292.

[7] Guéders, C., J. Van Roey, J. Gallant, and F. Coghe. “Simulation of Shock Wave Mitigation in Granular Materials by Pressure and Impulse Characterization.” In Proceedings of the 8th European LS-DYNA Users Conference, pp. 23–24. 2011.

[8] Zamani, Natasha, and Usama El Shamy. “Analysis of wave propagation in dry granular soils using DEM simulations.” Acta Geotechnica 6, No. 3 (2011): 167.

[9] Nguyen, Ngoc-Son, and Bernard Brogliato. “Shock dynamics in granular chains: numerical simulations and comparison with experimental tests.” Granular Matter 14, no. 3 (2012): 341–362.

[10] Neal, W. D., D. J. Chapman, and W. G. Proud. “Shock-wave stability in quasi-mono-disperse granular materials.” The European Physical Journal Applied Physics 57, no. 3 (2012): 31001.

[11] Nesterenko, Vitali. Dynamics of heterogeneous materials. Springer Science & Business Media, 2013.

[12] Rotariu, A-N., C. Dima, E. Trană, C. Enache, F. Timplaru, and L-C. Matache. “Uninstrumented Measurement Method for Granular Porous Media Blast Mitigation Assessment.” Experimental Techniques 40, no. 3 (2016): 993–1003.

[13] Kandan, K., Syed Nizamuddin Khaderi, H. N. G. Wadley, and V. S. Deshpande. “Surface instabilities in shock loaded granular media.” Journal of the Mechanics and Physics of Solids 109 (2017): 217–240.

[14] Liang, Tao, Xiao-Lu Zhu, Zhao-Cang Meng, Wen-Shan Duan, and Lei Yang. “Shock Waves in Nonlinear Granular Chain with Viscosity under Hertz Contact.” Journal of the Physical Society of Japan 88, no. 3 (2019): 034002.

[15] Rahmani, M., Oskouei, A. N., & Petrudi, A. M. (2020). Experimental and numerical study of the blast wave decrease using sandwich panel by granular materials core. Defence Technology.

[16] Cundall, Peter A., and Otto DL Strack. “A discrete numerical model for granular assemblies.” geotechnique 29, no. 1 (1979): 47–65.

[17] Ghaboussi, Jamshid, and Ricardo Barbosa. “Three-dimensional discrete element method for granular materials.” International Journal for Numerical and Analytical Methods in Geomechanics 14, no. 7 (1990): 451–472.

[18] Li, Xuefeng, Shibo Wang, Reza Malekian, Shangqing Hao, and Zhixiong Li. “Numerical simulation of rock breakage modes under confining pressures in deep mining: an experimental investigation.” IEEE Access 4 (2016): 5710–5720.

Biographies

Amin Moslemi Petrudi received his B.Sc. in Automobile Engineering from Allameh Amini university and M.Sc and PHD degrees in Mechanical Engineering from IHU University, Tehran, Iran. He has written numerous educational, research, and administrative records in his career. He is also capable and interested in research, design, modeling and simulation, impact and penetration mechanics, stress analysis.

Masoud Rahmani received his B.Sc. in Fluid mechanics from Malayer university and M.Sc and PHD degrees in Mechanical Engineering from IHU University, Tehran, Iran. He has written numerous educational, research, and administrative records in his career. He is also capable and interested in research, Simulation, Optimization and Nanofluid, stress analysis.

European Journal of Computational Mechanics, Vol. 29_4-6, 393–416.
doi: 10.13052/ejcm1779-7179.29464
© 2021 River Publishers