Numerical investigation of plane wave propagation in three soil models and a simple viscoelastic treatment of the boundary considering gravity using MPM

The importance of plane wave propagation in geotechnical engineering problems considering both small and large deformation is well understood but the theory of wave propagation in the latter is maybe not mature. In this paper, the mechanism of plane wave propagation in constitutive models such as linear elastic model, Mohr Coulomb model (MC) and elastoplastic density-dependent model (DDSM) considering finite deformations is extensively studied using the material point method. The classical laws of wave propagation and reflection in a continuum are found to provide reasonable explanations to the considered problems. The stress propagation mechanism in the three models considered is mainly related to the variation of the material wave impedance. An accurate silent boundary condition is found difficult even impossible to be determined for complicated models, which is verified by a simple deduction. Moreover, a simple viscoelastic silent boundary is proposed to consider the realistic effects of gravity in geotechnical engineering problems and its performance is compared with that of the classical viscous boundary for truncating a 1D linear elastic half-space.


Introduction
The study of wave propagation is involved in many practical problems in geotechnics. The theory of elastic wave propagation has already been developed as early as the nineteenth century, but stress wave in many practical problems is plasitc. Though limited analytical solutions have been proposed for plastic wave propagation problems, less work has been made to verify the numerical simulations with these analytical results [1]. Meanwhile, stress wave propagation for finite deformation problems have not gotten enough attention.
When the material domain of the problem is infinite, proper silent boundary conditions are necessary to be imposed on the truncated boundary which account for unobstructed propagation of waves from the near field to the far field, and minimize the reflection of waves on the computational domain boundary. Many formulations for silent boundaries have been proposed in connection with the finite element method. Despite their advantages, silent boundaries have limited applications in the material point method (MPM) [2].
MPM uses dual descriptions of Lagrange and Euler. A continuum is discretized into finite material particles that carry all material information. The motion and deformation states of the material are characterized by these material particles, while the Eulerian background grid is used to solve momentum equations and calculate spatial derivatives. During each time step, the material particles are solidly connected to the background grid. The material information is first mapped from material particles to grid nodes, and the momentum equations are solved on these grid nodes. Then, the information is mapped back to the material particles, and the material information on the grid nodes is set to be zero. By MPM, there are no convection terms, and the grid distortion is also avoided. Hence, it is suitable for dealing with large deformation problems.
In this paper, stress wave propagation mechanism for finite deformation problems in constitutive models such as linear elastic model, Mohr Coulomb model (MC) and Elastoplastic density-dependent model (DDSM) is investigated using numerical simulations. The difficulty in proposing an accurate silent boundary has been highlighted and verified by a simple deduction. As an alternate choice, the applicability of the traditional viscous silent boundary for finite deformation problems is investigated and a simple viscoelastic silent boundary is proposed considering the realistic effects of gravity. [3] based on researches by Harlow and Evans [4], Dawson [5], Brackbill and Ruppel [6], which combines the advantages of the Lagrangian and Eulerian methods. For the specific control equations, numerical processes and program implementations, readers are referred to the literature cited [7].

Material point method MPM was proposed by Sulsky et al
For standard MPM method, there is numerical noise caused by material particles crossing the background grid boundary. Zhang et al [8] proposed a new method for calculating the shape function gradient, by weighting the shape function gradient based on material particles in the standard material point method and the shape function gradient based on the background grid nodes in the fluid implicit particle method (FLIP) [6], to reduce the stress oscillation caused by material particles crossing the background grid boundary, which is called the dual domain material point method (DDMP).
In DDMP, the modified shape function gradient is as follows: (1) where the subscript I denotes the background grid node, p denotes the material particle, and j denotes the direction; is the shape function gradient of standard MPM, and is the shape function gradient based on the background grid nodes; is a weight function that is zero on background grid boundaries, to ensure the continuity of the modified gradient [9]. However, there is another problem for DDMP. The background grid node internal force is calculated by the modified shape function gradient, but the background grid node external force is calculated by the original shape function of standard MPM. This inconsistency will cause errors to some extent for the calculation of the deformation and stress near the material boundary. To overcome this small discrepancy, a modified shape function especially used to calculate the background grid node external force is proposed here: (3) (4) in which is the shape function gradient of standard MPM, and is the shape function gradient based on the background grid nodes.   For DDMP, the influence domains of a material particle and of a node are enlarged, which will induce the occurrence that a node near the material boundary will get node force from material particles but no mass. To overcome this problem, a simple but effective method is used to distribute the node force of a node with no mass to surrounding nodes with mass [8,10]. If this method is used, when information except node force is mapped from material particles to grid nodes, the original shape function of standard MPM is also effective, which means that only the calculation of grid node external force needs the modified shape function.

Silent boundary
The first implementation of silent boundaries in MPM was carried out by Shen and Chen [11]. They used the standard viscous boundary proposed by Lysmer and Kuhlemeyer [12] to simulate the film delamination. The viscous boundary is not only very easy to implement, but also quite accurate in many applications. Later on, several viscoelastic boundaries have also been proposed, but they are not the focus here.
The viscous damping stress on the boundary in 1-D case is as follows: (5) in which is the material density, is the velocity of the stress wave, and is the particle velocity. In MPM, the material particles carry all information such as density, stiffness, velocity, stress, etc., so as shown in figure 1, the boundary surface force is treated as the body force of the boundary layer, seen as the second term on the right side of the following equation: (6) in which is grid node external force; is the material particle mass; is body force; is the thickness of the boundary layer; is the specific boundary surface force (surface force over density); is the equivalent body force of the boundary surface force.
When applying viscous silent boundary in MPM, the outermost layer of material particles on the boundary is regarded as the boundary layer. The damping force of the silent boundary is approximately determined by the velocity, density and stiffness of the material particles in the boundary layer, then it is converted into the body force of the hypothetical boundary layer. This body force is then mapped to the background grid nodes and processed as node external force.
Equation (5) where the subscript 1 denotes the near field through which the stress wave passes first; 2 denotes far field through which the stress wave passes last. In this case, if we place a silent boundary just at the interface between these two materials, the accurate viscous boundary conditions should be equation (7) rather than 5. Equation (7) indicates that an accurate silent boundary should consider the material properties on both sides of the silent boundary. In more realistic geotechnical applications, the nolinearity, plasticity and large deformation all need to be considered, then an accurate silent boundary will be quite complicated, even impossible, due to the reflection of the stress wave caused by variation of material properties. As another choice, the variation of material properties near silent boundary is assumed to be small, and to investigate the applicability of the traditional silent boundary conditions in more realistic constitutive models will be meaningful.
In some dynamic cases, gravity may also have a certain significant impact. Here, a simple viscoelastic silent boundary is proposed considering the realistic effects of gravity: In the 1-D problem, in order to avoid tracking the distance from the silent boundary to the surface of soil, the above formula can be converted into: (9) in which is the initial material density; is the initial distance from the silent boundary to the surface of soil.

Numerical examples
The program of MPM is based on finite deformation theory, so all examples are actually finite deformation problems in this study. All numerical examples are 1-D here. Stress wave propagations in four material situations are considered. The first is in elastic material (EM). The second is in two different elastic materials, between which there is a interface (EM2M). MC model is used in the third one. An elastoplastic density-dependent model (DDSM) is used in the fourth one. In DDSM, the hardening of soil and the change of elastic modulus with variation of density can be considered, and the specific model is referred to the study of Zhang and Sun et al [9].
For convenience, some abbreviations are adopted. ID denotes that infinite material domain is adopted. SB denotes that the silent boundary condition shown as equation (5) is adopted. SB1 denotes that the silent boundary shown as equation (7) is adopted. SB2 denotes that the silent boundary shown as equation (9) is adopted. FB denotes that the fixed boundary is adopted. The depth of the 1-D bar are all 8 m for SB, SB1, SB2, FB, and the data in all figures are from points close to this depth or the boundaries.

Single elastic material
In order to approximate small deformation condition in MPM, the modulus is large and the load is relatively small. The elastic modulus is 5 Gpa, and the density is 2650 kg/m 3 . The input load is 100 kPa.
As shown in figure 2, the effect of SB is very close to that of ID, and the effect of SB is obviously better a lot than that of FB. The result of SB is a little different from that of ID, and this is maybe caused by the adoption of the finite deformation theory, which will be further verified in section 3.2.

Two different elastic materials
The second material is started from 8 m depth. The parameters for the elastic material through which the stress wave passes first are the same as section 3.1. For the second elastic material, the elastic modulus is 2.5 Gpa, and the density is 1325 kg/m 3 . The input load is 100 kPa.
Calculated by equation (7), the stress on the interface should be 67 kPa in small deformation condition. As shown in figure 3, the results of SB1 and ID are close to this value, while that of SB deviate form this value a lot. It means that, if the material properties change along the stress wave propagation path, the accurate silent boundary condition will be very complicated. It is difficult even may be impossible to get a totally accurate silent boundary condition in more realistic soil models, because the material properties on both sides of the silent boundary should be considered.
According to equation (7), if of the second material is smaller than that of the first material, then the total stress on the interface will be smaller than the incident stress because the superposition of the reflection stress to the incident stress. If besides the incident stress the refection stress is absorbed by the silent boundary at the same time, then the total stress close to the silent boundary will increase.
is called wave impedance. So the result of SB is larger than that of ID in figure 2, and the result of SB1 is larger than that of ID in figure 3, due to that the finite deformation theory is adopted and the wave impedence in front of the stress wave is smaller than that behind the stress wave.
1-D elastic stress wave theory can be used to discuss the error of the silent boundary due to the absorption of the reflection stress wave. This error is: (10) (11) in which is the incident stress; is the initial stress of one point. If the variation of the material properties and the difference between the incident stress and the initial stress are all very small, then this error will be very small. This law should be the same for more complicated models. If the energy is absorbed by the plastic deformation, and the silent boundary is placed enough far from the load end, then and will be very small, and the silent boundary will be more accurate.

MC model
The Young's modulus is 50 MPa; the density is 2650 kg/m 3 ; the Poisson's ratio is 0.3; the friction angle is 30°; the other parameters are zero. The input load is 100 kPa. As shown in figure 4, the result of SB is also close to that of ID, and obviously better a lot than that  It should be noted that, the input loads for section 3.1, 3.2, 3.3 are all relatively small, so the deformations are close to small deformation, and the variation of the material properties may be not so obvious. Due to space limitations, the influence of large variation of the material properties for single elastic material and MC model is not investigated here and that will be shown directly in simulation of DDSM model.

DDSM model
The specific model is referred to the study of Zhang and Sun et al [9]. The parameters are shown in table 1, in which is the void ratio at pressure 0, is called grain hardness and is a parameter reflecting the grain skeleton characteristics of soil, m is called stiffness pressure correlation index related to the properties of soil, is the ratio of the unloading modulus to the loading modulus, is Poisson's ratio, c is cohesion, is friction angle, and is dilatancy angle. Because energy is absorbed very quickly in this model, the input load is 40 MPa. As shown in figure 5, the stress propagation is quite complicated. But the peak of the stress and the peak of the plastic strain energy occur at the same time, it means that the stress after the peak will not produce a substantial impact, so we can focus on the peak of the stress only.
As shown in figure 6, though the result of SB deviate from that of ID to some extent, but the result of SB still quite better than that of FB.

Gravity
The single elastic material is adopted and the parameters are the same as that in section 3.1. As shown if figure 7, if gravity is considered, the result of SB will occur rigid motion, while SB2 proposed in this study behave well.

Conclusions
In this paper, the mechanism of plane wave propagation in constitutive models such as linear elastic model, MC model and DDSM model considering finite deformations is extensively studied using MPM. The stress propagation mechanism can be regarded as mainly related to the variation of the material wave impedance. A totally accurate silent boundary condition is found difficult even impossible to be determined in complicated models, which has been verified by a simple deduction. The error of the traditional silent boundary in models such as MC and DDSM, due to the absorption of the refection stress caused by variation of the material properties, will be very small if the silent boundary is placed enough far form the load end (of course not too far). A simple viscoelastic silent boundary is proposed considering the realistic effects of gravity in geotechnical engineering problems and it performs well. The results presented in this work are a result of preliminary investigations only and detailed work will follow later.