Abstract
From the microperspective, this paper presents a model based on a new type of noncontinuous theoretical mechanical method, molecular dynamics (MD), to simulate the typical soil granular flow. The Hertzian friction formula and viscous damping force are introduced in the MD governing equations to model the granular flow. To show the validity of the proposed approach, a benchmark problem of 2D viscous material flow is simulated. The calculated final flow runout distance of the viscous material agrees well with the result of constrained interpolated profile (CIP) method as reported in the literature. Numerical modeling of the propagation of the collapse of three-dimensional axisymmetric sand columns is performed by the application of MD models. Comparison of the MD computational runout distance and the obtained distance by experiment shows a high degree of similarity. This indicates that the proposed MD model can accurately represent the evolution of the granular flow. The model developed may thus find applications in various problems involving dense granular flow and large deformations, such as landslides and debris flow. It provides a means for predicting fluidization characteristics of soil large deformation flow disasters and for identification and design of appropriate protective measures.
1. Introduction
As a kind of typical granular material, soil causes the common flow forms mainly for earthwork excavation, flow of soil slope by filling, and sudden slip on a weak foundation soil. The large deformation disasters caused by the flow of soil particles have received considerable attention and researches both at home and abroad.
Molecular dynamics originated in the 50s and began to be widely concerned in the mid-70s. In 1957, the state equations of gas and liquid were studied, firstly using molecular dynamics in the hard sphere model, thus setting a precedent for studying the macroscopic properties in the molecular dynamics method [1]. Molecular dynamics started to develop after 70s, gradually perfected its system in the late 80s and got rapid development along with the flourish of computer hardware in the 90s. Molecular dynamics is a connecting link between macroscopic characteristics and microstructure, which can make certain microscopic interpretation for theoretical analysis and experimental researchdifficult to understand. Especially in chemistry, biology, physics, and material science or related fields, it gets better application [2–4]. In recent years, molecular dynamics is also applied to the research on granular flow of geomaterials. The inducing factors of shallow landslide were studied using the molecular dynamics method [5]. Factors affecting the occurrence of microsliding in the sand pile were discussed with the molecular dynamics method [6]. Molecular dynamics simulations are particularly well suited to investigate the effects of materials properties of liquid-solid interfaces on flow boundary conditions [7].
Molecular dynamics method, which is a combination of physics, mathematics, and chemistry, mainly relies on Newtonian mechanics to simulate the movement of molecular system and can simulate the movement of a single particle in a large number of particle collection system (solid, gas, or liquid). Its key is the concept of movement, that is, to compute the evolution of the position, velocity, and orientation of particles over time. The particle in molecular dynamics can be atom, molecule, or larger set of particles [8]. The basic principle of molecular dynamics is the application of Newton’s classical mechanics to firstly calculate the trajectories of all particles in the physical system and then by using relevant statistical methods to calculate mechanics, thermodynamics, and dynamics characteristics of the system. That is to say, firstly, the system of particles is abstracted to many interacting particles, and the interaction force between particles is obtained by using the potential energy function derivation method. And then, the acceleration and velocity of each particle are calculated according to the Newtonian mechanics, and the particles have the corresponding trajectories after a certain integral steps, which must be preserved by setting a certain time interval. Finally, calculation results of the interested quantities are extracted based on statistical principle [9–11].
At present, the molecular dynamics simulation researches mostly focus on physics, biology, and chemistry, but the application in the field of geotechnical engineering is rare. The molecular dynamics approach is suitable to model granular material and to observe the trajectory of a single particle, so as to possibly identify its dynamical properties. Thus, the molecular dynamics method is tried to be applied to the granular flow characteristics research, and a 3D MD model that simulates granular flow is presented in this work, which is dedicated to provide a strong scientific basis for solving geotechnical engineering problem.
2. Numerical Approach
The accuracy of the model depends on the accuracy of description of external forces acting on the particle and the treatment of collisions. The particles can be either rigid and/or soft. A soft particle model [12] is used in analyzing particle-particle and particle-wall interactions in this paper. The particles are considered to be spherical and identical. The collisions are assumed to be central and involve both linear elastic and damping forces. In the study, 3D ordered granular packings made up of spherical particles are considered. The motion equation for the sphere is written in a general form aswhere is the displacement, is the force, is the mass, and is the total number of particles. is the composition of contact force between spheres j and i, where the interaction force between any contacting spheres under compression is given by the Hertzian law [13]:where is the overlap distance of two particles and the and stand for the radius of the particles i and j, respectively. The expression of is defined as [13]where and are the elastic constants for normal contact and tangential contact, respectively; and are the viscoelastic damping constants for normal contact and tangential contact, respectively; is the effective mass of two particles of mass and ; is the tangential displacement vector between two spherical particles which is truncated to satisfy a frictional yield criterion; is the unit vector along the line connecting the centers of the two particles; and and are the normal and tangential components of the relative velocity of the two particles, respectively.
Viscous damping force calculation formula of the particle is based on the Stokes drag law of fluid mechanics [14], which can be expressed aswhere is the dynamic viscosity of the frictional fluid, is the diameter of particle, and is the velocity.
At present, in the field of molecular dynamics, the integral of the motion equation mainly adopts Verlet algorithm [15]. The integral algorithm putted forward by Verlet is the most widely used in the molecular dynamics. The numerical simulations are performed using the molecular dynamics package LAMMPS [10] in the work, by considering spheres as point masses connected by nonlinear springs.
The calculation formulas of Verlet algorithm are as follows:
3. Validation of the MD Model
The performance of the new model will be verified in this section. It is essential to create an input script containing the desired commands before running LAMMPS. The Verlet algorithm is applied to the integral of motion equation in the input scripts. When numerical calculation starts, LAMMPS reads the input script, and the displacement, velocity, and acceleration of every granular particle are calculated and updated, thus obtaining the evolution trajectory of the whole calculation system over time. A two-dimensional simulation of viscous material flow was conducted to verify the reliability of the MD method.
Viscous materials flow on the plane under its own gravity, and a literature [16] has detailed calculation and analysis on this flow process using constrained interpolated profile (CIP). CIP method was developed for prediction of large deformations associated with a geomaterial flow, so it can be used as a contrast with the MD method in this work. The model [17] is shown in Figure 1.
Based on the molecular dynamics model mentioned above, the viscous material flow process was simulated. In MD simulation, the calculation area of 50 m long and 20 m wide and the viscous material size of 12 m long and 10 m wide were set up. The particles were idealized as circle with ignoring particle size distribution. The boundary conditions were fixed, and the interaction between particles in the system was calculated according to the Hertzian friction formula. The initial velocity was set to zero.
Calculation parameters and simulation results are shown in Table 1 and Figure 2, respectively. The molecular dynamics simulation results are compared with the CIP results of the existing literature [16], to verify the validity and correctness of molecular dynamics algorithm in particle flow simulation. Contrast can be seen in Figure 3, and the red parts are the MD numerical simulations of flow profile, while black lines show CIP method calculation configurations in the literature [16].
(a)
(b)
(c)
(a)
(b)
(c)
Figure 3 shows that the MD results are consistent with the CIP results. The literature [16] discussed in detail the applicability and accuracy of the CIP method in describing the large deformation and fluidization motion characteristics of the viscous material. Therefore, the applicability and accuracy of the MD program in the description of fluidization motion characteristics of the viscous material have also been verified and validated.
4. Application of MD Modeling to Collapse Test and Numerical Simulation for Sand Column Collapse
This section discusses the three-dimensional simulation of the propagation of sand flow as the application of MD modeling to real large deformation of geomaterials. The motions of sands after collapse were represented in the MD simulation. This study can produce preliminary results for the prediction fluidization characteristics of large deformation flow disasters of geomaterials.
4.1. Sand Column Collapse Experiment
The collapse test of three-dimensional axisymmetric sand columns was carried out on the basis of the apparatus developed independently, which was made of steel with the inner diameter of 10 cm, the height of 10 cm, and the thickness of 1 cm.
Firstly, the collapse device was placed on a smooth horizontal plane to ensure that there was no gap between the experimental device and the contact surface. Secondly, the dry sand was slowly sprinkled into apparatus to keep uniform distribution. When the sand sample reached the height of the design of 10 cm, it was stopped and made to stand for a period of time, until the sand was compacted under its own gravity. Then, the camera was opened, and the collapse device was vertically and quickly removed. Finally, configuration images of sand were collected and analyzed. In the experiment, the frictional effects between individual grains on granular flow play a role only in the last instant of the flow, as it comes to an abrupt halt.
The final configuration of sand after collapse is conical in the experiment, as shown in Figure 4. The average value of the cone diameter is 28.87 cm. The time-history curve of the cone diameter in the slump test was obtained and shown in Figure 5. It shows that the cone bottom diameter has a trend of increase with time. The slope of the polyline indicates the speed of growth and the diameter growth rates are obviously different in the whole collapse process. As the evacuation of the test device has a certain friction effect on sands in the experiment, the forces of sands are not stable in the whole process, which may bring certain influence on the flow rate and the final configuration of the sands.
4.2. Simulation of Granular Flows during the Collapse of Sand Columns
The developed and validated MD model was used to simulate the collapse propagation of an initially vertical 3D axisymmetric sand column in the work. The MD input file for 3D simulation of sand columns collapse was compiled, and then the diameter of the bottom circle formed by the final flowing configuration was calculated using the LAMMPS package. In the proposed model, the sand materials were discretized into a series of MD particles with a certain diameter. In the numerical simulation, sands were regarded as the uniformly distributed spherical particles, the calculation area was a cylinder with the diameter of 10 cm and the height of 10 cm, and the initial velocity of every sand particle was set to zero. The parameters used in the MD simulation were listed in Table 2. Calculating diameter was set 3 times of the expectation value of the particle diameter about 2 mm. The value of elastic constant for normal contact of the sand was 10 MPa, and the value of elastic constant for tangential contact was set two-sevenths of the former. The damping constant for normal contact of the sand was set as 0.1, and the damping constant for tangential contact was half of the former. The time step was 1.0 × 10−5 s, and there were about 118,700 MD particles used to represent the sand column in the simulations conducted. The final runout distance was investigated and compared with the experimental observation, and the result was presented subsequently.
Comparisons of MD numerical simulation results and the experimental results are shown in Figure 6. Excellent agreement is observed between the experiments and 3D MD simulations conducted with respect to the flow distance. Obviously, Figure 6 illustrates that there is a small error between the MD simulation results and the experimental results. The possible reason is that the MD numerical calculation process has ignored the friction effect on sand of experimental device.
Thus, the MD method can describe well the granular flow of sand columns after collapse, and the correctness and applicability of the model and algorithm of molecular dynamics are verified again.
5. Conclusions
Application of the MD method to the simulation of soil granular materials under large deformation is the subject of this paper. The Hertzian friction formula and viscous damping force are implemented in the MD framework to model the granular flow.
The model developed is validated by the collapse of viscous material flow as reported in the available literature. Excellent agreement is observed between the MD model simulations and CIP simulations with respect to the flow distance, thus verifying the performance of the MD modeling scheme.
Simulation of the collapse of 3D axisymmetric sand columns is also conducted. Numerical results for the granular flow pattern and final runout distance are in good agreement with the experimental observations. The flow velocity of the simulated deformed region within a sand column during the collapse is bigger than the experimental observation. The likely cause is that MD numerical simulations ignore the influence of the evacuation of collapse device on the sand, whereas the collapse device has created the friction on the sand particles when it is removed in the experiment process.
This study indicates that the MD model developed can be used to effectively simulate large deformation of granular materials, and geomaterials in general, if proper calculation models are implemented. This is due to the fact that the MD method is a particle-based mesh-free numerical method. In the MD method, the variables of particles are calculated through the integration procedure. If an appropriate integration step is chosen, the MD model developed is capable of presenting many discrete characteristics of dense granular flow for granular materials under large deformation. Thus, the MD model may find extensive applications in various problems involving granular flow such as landslides and debris flow.
Conflicts of Interest
The authors declare that they do not have any commercial or associative interest that represents a conflicts of interest in connection with the work submitted.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (NSFC) (Grant no. 51608316), the Research Fund of the Key Laboratory of Transportation Tunnel Engineering, Ministry of Education (Grant no. TTE2014-05), and the Traffic Science and Technology Program in Shanxi Province (Grant no. 2017-1-20). The authors would like to express their gratitude for the financial support.