Volume 8, Issue 3

September 2023

Special Contribution

Numerical Simulation of Penetration into Ground
Yosuke Higo
Professor, Graduate School of Management (double appointment)
Graduate School of Engineering, Dept. Urban Management
Kyoto University

1. Introduction
When structures such as sheet piles and piles are driven into the ground, the ground around the structures is significantly
disturbed by large shear forces due to friction between the ground and the structures. In the case of pile driving, large
deformation also occurs at the pile tip. It is very important to accurately analyze such large deformation of the groundand to evaluate the impact of structure penetration on the surrounding ground. This article reviews the Material PointMethod (MPM), which is one of the mesh-free deformation simulation methods for soils. MPM has been rapidly
developed as a method for simulating large deformation of solids and is spreading internationally. Applicability of MPM
to numerical simulation of penetration into ground is shown through numerical examples of sheet pile penetration andfooting penetration simulations.

2. Material Point Method
Numerical simulation techniques within a continuum mechanics framework are classified into three primary groups
based on their definition of the element control domain. The initial category is the mesh-based approach, represented
by the Finite Element Method (FEM). In FEM, the element domain is determined by nodes and their connectivity,
distinctly outlining the control domain through node coordinates. Though adept for boundary value problems, FEM canface mesh-tangling during severe deformation due to the connectivity of elements. The second category involves particlebased method, where material properties gather around representative points. Interpolation functions define the circular(2D) or spherical (3D) control domain shape. Well-known within particle-based methods are the Moving Particle Semiimplicit (MPS) method (e.g., Koshizuka and Oka, 1996) and Smooth Particle Hydro-dynamics (SPH) method (e.g.,Monaghan and Lattanzio, 1985), which describe the material domain through particle assemblies. Unlike FEM, particlebased methods lack clear material boundaries, necessitating special treatment of boundaries. Conversely, particle-based
methods can avoid mesh-tangling since individual particle domains, formed by interpolation, remain disconnected. The
third type is the grid-particle hybrid method, typified by the Particle-in-Cell (PIC) method (Harlow, 1956), Fluid-ImplicitParticle (FLIP) method (Brackbill and Ruppel, 1986) for fluid mechanics, and Material Point Method (MPM) (Sulsky et al.,1994) for solid mechanics. These methods also depict the material domain through particle assemblies. While similar toother particle-based methods in simulating objects through arranged particles, grid-particle hybrids offer clear boundary
condition definition by applying conditions to grids. The boundary setup in hybrid methods resembles FEM, facilitatingeasier configuration compared to particle-based approaches. In these methods, material properties center on materialpoints, while governing equations are solved at grid nodes. Interpolation functions aid in exchanging physical quantities
between material and grid nodes.
Sulsky et al. (1994) initially introduced MPM as a derivative of the PIC method for solid mechanics, termed the "originalMPM" in this context to differentiate from other MPM derivatives. Ongoing MPM evolution mainly centers on integrating
control domain deformation. A novel method applied to geomaterial deformation, termed Arbitrary Particle Domain
Interpolation (APDI) permits flexible material point control domain definitions (Kiriyama and Higo 2020). Fig. 1 illustrates
the relationships among various 2D and 3D MPM implementations, defining their control domains. The original MPM
represents a material point as a simple point (Fig. 1(a)), but this leads to numerical oscillations when material points crossgrid boundaries. To mitigate this, Bardenhagen and Kober (2004) proposed the Generalized Interpolation Material Point(GIMP) method. GIMP manifests in uniform/unchanged (uGIMP) and contiguous particle (cpGIMP) forms. In cpGIMP (Fig.
1(c)), the control domain updates with deformation, while Sadeghirad et al. (2011, 2013) further advanced GIMP into theConvected Particle Domain Interpolation (CPDI) method, accommodating shear and rotational control domain
deformations. CPDI divides into CPDI1 (Fig. 1(d)), with parallelogram deformations, and CPDI2 (Fig. 1(e)), withquadrilateral deformations. Extending to three dimensions involves expanding from two dimensions for the original MPM,GIMP, and CPD1 methods. CPDI2's three-dimensional formulation necessitates a distinct approach, using direct
integration for volume integration. The APDI method circumvents this complexity through numerical integration, offering
a comprehensive formulation for both 2D and 3D, along with adaptable control domain deformations. CPDI2 uses3

Volume 8, Issue 3

September 2023

quadrilateral control domain deformation, whereas APDI accommodates diverse control domain deformations.
The governing equations of MPM are fundamental laws of mechanics: the conservation of mass and the linearmomentum. The mass conservation always holds because the material point possesses mass throughout the calculations.
Hence the conservation of the linear momentum, i.e., the equation of motion, is only one equation to be solved. Notethat governing equations for liquid phases and interaction with the solid phase, e.g., Darcy’s law, are needed whensimulating hydro-mechanical coupled behavior of soils as porous media. The algorithm of MPM is schematically shownin Fig. 2. The continuum body is divided into multiple subdomains and represented by material points with mass of thesubdomain. The history-dependent variables such as stress and strain are stored and delivered by material points. Hence,
history-dependent constitutive equations can be introduced. The governing equations are solved on the background
Eulerian computational grid and thus MPM does not face to mesh tangling problem even when the severely large
deformation. Similar to FEM, the boundary conditions are imposed on the computational grids. Although the grids
containing the particles can roughly define the boundaries of a continuum body, exact representation of the boundaries
need coupling with element-based method.

2DPointSquare

Rectangular

Parallelogram

Quadrilateral(CPDI2)

TriangleHexahedron

Tetrahedron

3DPointCube

Rectangular

Parallelopiped(a) Original MPM(b) uGIMP(c) cpGIMP(d) CPDI1(e) APDI

Fig. 1. Variations of particle domain for Material Point Method

Fig. 2. Numerical algorithm of Material Point Method

3. Steel Sheet Pile Penetration
Numerical simulation of steel sheet pile penetration was conducted using a displacement-control method for simulating
quasi-static press-in penetration and a stress control method for simulating dynamic hammer-blow penetration. Theparticle configuration and the boundary conditions are shown in Fig. 3. The thickness of the steel sheet pile is describedby four particles allocated in two grids, which is much thicker than the real sheet piles, to avoid numerical instability whenusing smaller number of particles. The larger number of particles in one computational grid generally makes the accuracyof MPM simulation higher. The steel sheet pile is assumed to be a linear elastic body whose Lame constants 𝜆 and 𝜇
are 4.2×107 kN/m2 are 2.80×107 kN/m2, respectively. Although the second moment of area is larger than that of the real
sheet pile, the set of elastic constants is determined to satisfy the bending stiffness of a real steel sheet pile. Consequently,
the axial stiffness is underestimated. The model ground is an elasto-plastic body described by Mohr-Coulomb type yieldfunction and Drucker-Pager type potential function. The internal friction angle and the cohesion are 30 degree and 8.5
kN/m2, respectively. Particles of the ground in the vicinity of the sheet pile are modeled as interface particles whose
internal friction angle and cohesion are assumed to be 20 degree and zero, respectively.
4

Volume 8, Issue 3

September 2023

The loading conditions are shown in Fig. 3. In the displacement-controlled method, a constant displacement of 0.02
m/s was applied to the steel particles above the ground surface to penetrate the sheet pile. In the stress-control method,
the top four steel particles were subjected to linearly varying surface forces in the form of triangular waves at intervalsof 0.04 s to a maximum of 100 kN per unit depth, for a total of 400 kN, and the same surface forces were applied again
after 0.56 s, repeated 19 times. The results of the analysis for a single blow were conducted beforehand, and the blow
interval was determined as the time required for the ground vibration caused by the blows to settle down. The numberof blows was set to be sufficient for the steel head to penetrate below the ground surface.
(m)

0.02 m/s

Static penetration(displacement control)

(sec)

(kN)
19 blows

Dynamic penetration(force control)

………
(sec)

Fig. 3. 2D sheet pile penetration model (unit: m) and the loading conditions.

(a) Static press-in penetration(b) Dynamic hammer-blow penetration
Fig. 4. Numerical simulation of sheet pile penetration into ground. Maximum shear strain distribution for the static anddynamic penetration cases.

(a) x-direction(b) y-direction
Fig.5. Distribution of acceleration for the dynamic penetration case just after the first blow (0.12 second).

5Volume 8, Issue 3

September 2023

Fig. 4 show the maximum shear strain distribution at the end of penetration for the static and dynamic penetration. Theresults show that both the displacement-controlled and stress-controlled methods can simulate the behavior of steel
penetration into the ground with large shear strain around the sheet pile. In the case of static penetration, the sheet piledrags the surrounding soil into the ground as the sheet pile penetrates, resulting in a large settlement. On the other hand,in the case of dynamic penetration, the steel penetration is observed without significant drag of the ground as observedin the case of the static press-in penetration.
The acceleration distributions in the x and y directions for the first blow in the dynamic penetration analysis are shownin Fig. 5. In dynamic penetration, the acceleration is distributed alternately positive and negative, and the acceleration
propagates in the direction away from the sheet pile. In other words, the impact by the blow generates sparse and dense
waves in the ground. While the ground is sparse near the steel, the confining pressure on the steel is reduced, which alsoreduces the frictional resistance between the ground and the sheet pile, which suppresses the dragging of thesurrounding ground. On the other hand, in static penetration, the steel is constantly confined by the surrounding ground
during penetration. Thus the frictional resistance is not reduced and the shear surface due to the frictional force reaches
the ground surface as shown in Fig. 5. This causes significant drag of the surrounding ground over a wide area.
To evaluate even the stress level of sheet pile, much more particles must be placed in the thickness direction of the sheet
pile, which is computationally too expensive. Although such an analysis is possible with a large-scale analysis using parallel
computation, it will be difficult to implement in practice in the near future. However, it is possible to evaluate the behaviorof the surrounding ground at least in practice by correctly evaluating the interaction between the ground and thestructure.

4. Footing Penetration
Computational cost of three-dimensional numerical simulation by an explicit numerical integration scheme is too
expensive to be done. A practical implementation scheme, the particle-element coupled method has been proposed. In
this scheme, with the GIMP interpolation function, highly deformable regions are represented by material points. With
the APDI interpolation function, material points are used, as if they were FEM elements, in the rest of the highly
deformable regions. This combination of the GIMP and APDI methods makes it possible to divide the material domain
into spatial discretization with different sizes. The particle-element is able to simulate extreme deformation by GIMP
particles as well as to evaluate boundary values by APDI elements for the area of focus, which is sometimes difficult inpurely particle-based methods due to the lack of an element boundary. The boundary between the particle and theelement domains does not require any special treatments, but simply use interpolation function which satisfies partitionof unity for each domain.
The three-dimensional model of the shallow foundation bearing capacity problem is shown in Fig. 6. The space aroundthe footing and down to the bottom of the base ground is modeled with GIMP particles in a cubic manner. Each
computational grid is a cube, 0.2 m on each side, and 8 particles are allocated per cell for the GIMP particles. Considering
the symmetry of the scenario, a quarter of the domain is modeled. The normal direction of the lateral boundary is fixed(vertical roller condition), and the bottom surface is fixed. The other conditions for the three-dimensional analysis arelisted in Table 1. Vertical load is incrementally applied to the base ground from 100 kPa to 600 kPa.
The snapshots of maximum shear strain distribution are shown in Fig. 7. At 100 kPa, the maximum shear strain is
distributed like a stress bulb in the Boussinesq solution based on the elastic theory. At 200 kPa, the maximum shear strain
is concentrated at the outer edge of the footing, forming a punching wedge beneath the footing. This strain concentration
at the outer edge of the footing is identical to the theoretical solution for cohesive material. In the regime beyond 300
kPa, the base ground reaches the slip failure mode and the resisting slip line spreads in space. Settlement then increases
continually to 400 kPa and surface deformation becomes visible in Fig. 7(d). However, there is no visual discontinuity inthe base ground at this stage. The loading block continues to penetrate into the base ground up to 500 kPa and 600 kPa,
where discontinuity between the loading block and the base ground becomes visible (Fig. 7(e) and (f)).
These results demonstrate that the particle-element coupled method is applicable to three-dimensional problems,
reducing the implemental complexity and computational cost by dividing the domain into two parts (a large deformation
domain and a small deformation domain). The computational cost is reduced by implementing the GIMP particles in thelarge deformed domain and the APDI particles in the small deformation domain. Use of a full particle model for the threedimensional problem would result in an extremely high computational cost as compared with the two-dimensional model,
but the particle-element coupled method makes it possible to perform the three-dimensional simulation at a reasonable
computational cost.
6Volume 8, Issue 3

APDI
Loading Block48
48
4uGIMP

Items
ValuesParticles per cell
8Dimensions(H x W x D)
6 m x 12 m x 12 m
Width of cell
0.2 mTime increment
0.00025
Damping factor
0.8
Convergence criteria
0.0001
Number of uGIMP interpolation Material Points
Nearby Ground64,000Number of APDI interpolation Material Points
Loading Block
8,000
Surrounding Ground
24,580Fig. 6. 3D element-particle coupled footing penetration model(unit: m)

00.010

September 2023

Table 1. Simulation conditions

00.047

Elastic bulb

Wedge formation(a) 100 (kN/m2)0(b) Pressure=200 (kN/m2)

0.15800.838

Vertical concentration(c) Pressure=300 (kN/m2)0(d) Pressure=400 (kN/m2)

3.11

0

5.47

Penetration(e) Pressure=500 (kN/m2)(f) Pressure=600 (kN/m2)

Fig. 7. 3D numerical simulation of footing penetration into ground. Maximum shear strain distribution at increasing
7

Volume 8, Issue 3

5. Summary and Future Perspectives

September 2023

loading pressures.

The numerical simulations of penetration into ground successfully performed by MPM are demonstrated. To evaluate
even the stress level of steel, multiple particles must be placed in the thickness direction of the steel, which requires avery large number of particles and is computationally too expensive. Although such an analysis is possible with a largescale analysis using parallel computation, it will be difficult to implement in practice in the near future. However, it ispossible to evaluate the behavior of the surrounding ground at least in practice by correctly evaluating the interaction
between the ground and the structure.
In general, design and analysis are performed assuming ideal conditions after sheet pile and pile driving. In addition toMPM, several other methods have been developed that can be applied to large deformation simulations. Among them,
MPM is similar to FEM in that the governing equations are solved on a background grid, and users of FEM, which iscurrently the most popular numerical analysis method, can use MPM relatively easily.

Acknowledgements
The author expresses sincere gratitude to the following three colleagues who performed the numerical simulations
shown in this article: Dr. Takatoshi Kiriyama (Institute of Technology, Shimizu Corporation), Mr. Yudai Takegawa (Railway
Technical Research Institute), and Gengo Kamamori (mater course student, Kyoto University).

References
Bardenhagen, S.G., Kober, E.M., The generalized interpolation material point method. Comput. Model. Eng. Sci. 5 (6),447–495, 2004.
Brackbill, J.U. and Ruppel, H.M., FLIP: A method for adaptively zoned, particle-in-cell calculations in two dimensions,
Journal of Computational Physics, 65, 2, 314-343, 1986.
Harlow, F.H., A machine calculation method for hydro-dynamic problems, Los Alamos Scientific Laboratory report LAMS,
1956.
Kamamori, G., Higo, Y., Takegawa, Y., Miyoshi, T., Analysis of subsoil behavior during static and dynamic steel penetrationby Material Point Method, Proc. The 57th JGS annual conference, 21-6-1-05, 2022. (in Japanese)
Kiriyama, T. and Higo, Y., Arbitrary particle domain interpolation method and application to problems of geomaterialdeformation, Soils and Foundations, 60 (6), 1422-1439, 2020.
Koshizuka, S. and Oka, Y., Moving particle semi-implicit method for fragmentation of incompressible fluid, Nuclear Scienceand Engineering, 123, 421–434, 1996.
Monaghan, J.J. and Lattanzio, J.C., A refined particle method for astrophysical problems, Astronomy and Astrophysics,
149, 1, 135-143, 1985.
Sadeghirad, A., Brannon, R.M., Burghardt, J., A convected particle domain interpolation technique to extend applicabilityof the material point method for problems involving massive deformations, Int. J. Numer. Meth. Eng. 86 (12), 1435–1456, 2011.
Sadeghirad, A., Brannon, R.M., Guilkey, J.E., Second-order convected particle domain interpolation (CPDI2) withenrichment for weak discontinuities at material interfaces, Int. J. Numer. Meth. Eng. 95 (12), 928–952, 2013.
Sulsky, D., Chen, Z., Schreyer, H.L., A particle method for history-dependent materials, Comput. Methods Appl. Mech.
Eng., 118, 179–196, 1994.

8