Impact of Pile Punching on Adjacent Piles: Insights from a 3D Coupled SPH-FEM Analysis

: Pile punching (or driving) affects the surrounding area where piles and adjacent piles can be displaced out of their original positions, due to horizontal loads, thereby leading to hazardous outcomes. This paper presents a three-dimensional (3D) coupled Smoothed Particle Hydrodynamics and Finite Element Method (SPH-FEM) model, which was established to investigate pile punching and its impact on adjacent piles subjected to lateral loads. This approach handles the large distortions by avoiding mesh tangling and remeshing, contributing greatly high computational efficiency. The SPH-FEM model was validated against field measurements. The results of this study indicated that the soil type in which piles were embedded affected the interaction between piles during the pile punching. A comprehensive parametric study was carried out to evaluate the impact of soil properties on the displacement of piles due to the punching of an adjacent pile. It was found that the interaction between piles was comparatively weak when the piles were driven in stiff clays; while the pile-soil interactions were much more significant in sandy soils and soft clays. when the clear spacing is 3 times the pile diameter. Tilt of the adjacent pile was calculated as the ratio of head displacement relative to the tip displacement of the pile to pile length. In all the cases, it was observed that the tilt of the adjacent pile is insignificant. The maximum tilt of 0.00304 (i.e., about 1/329) was obtained for the very loose sand which has elastic modulus of 10 MPa. The results show that tilt of pile decreases when the pile is installed in a dense or hard soil.


Introduction
Piles are commonly used as foundations for many major structures in civil engineering to transfer the heavy loads, for which shallow foundations may not be economical and feasible. In engineering practice, piles are often designed only for carrying vertical loading, as typically the vertical loads are significantly larger than the horizontal loads such as wind loading. However, piles can also be subjected to the lateral loads from the surrounding soil due to construction activities [1,2]. The lateral soil movement generated by the installation of a pile close to exiting piles will induce additional deflection and bending moment to the adjacent existing pile. Thus, the lateral response of pile foundation is also essential in the designing of structures where lateral dynamic loads exist. Unlike the axial load capacity of a pile, the determination of its lateral load capacity is much more complicated because the soil-pile interaction affects the pile deflection [3,4].
When the punching hammer strikes a pile head, a stress wave is generated within the pile that travels along the pile, during which a part of the energy is transmitted into the soil at the soil-pile interface [5][6][7][8]. Thus, the pile punching affects the surrounding area where piles are installed, and the adjacent piles can be displaced from their original positions. Additionally, the vibrations induced by pile punching can damage structures and cause discomfort to the people in the proximity of pile punching. Thus, the prediction of ground vibration from pile punching, and the study of its impact on adjacent piles, are crucial in preventing potential damage on the nearby environment and structures.
Previous studies have focused on the driving efficiency of piles, and only few investigations of vibrations due to pile punching, and their effects on the nearby structures, are available [9][10][11][12][13]. Nowadays, Finite Element Method (FEM) analysis has become a promising approach for studying the problems in soil-pile interaction. It is known that the soil in the vicinity of the pile can be subjected to large deformations (as a result of pile penetration), and FEM analysis, therefore, should have the ability to consider large deformations [14]. Thus, researchers have used different techniques to simulate the pile driving, such as lumped parameter models, Material Point Method (MPM), and continuum FEM models using Arbitrary Lagrangian-Eulerian (ALE) method. However, in most previous FEM analysis, the installation of piles has not been explicitly modeled and two-dimensional (2D) axisymmetric models are often used [8,[15][16][17][18][19].
Apparently, these 2D analyses could not incorporate the radial and three-dimensional componentinteractions. Thus, they are not well-suitedfor understanding the pile-soil-pile interaction in a real environment. As such, a 3D FEM analysis is needed to unveil the real interaction mechanism. However, these 3D models have rarely been available in literature, given the 3D FEM analysis requires a considerable computational effort for generating input and interpretation of results. Also, the impact mechanism underlying the dynamic interaction between adjacent piles in the process of pile punching is still not clear, although some investigations are available [20,21].
Another limitation of FEM in the application of large deformation problems is that the use of conventional Lagrangian meshes will result in mesh tangling, leading to severe numerical instabilities. Smoothed Particle Hydrodynamics (SPH) method has a strong ability to solve dynamic problems involving large deformation. On the other hand, it is not as good as the FEM in terms of computational time and boundary conditions. In this regard, coupled SPH-FEM method can be effectively used of two kinds of algorithms for the simulation of large deformation problems by eliminating the limitations in those two algorithms. Today, vast number of FE codes are available that are capable of analyzing challenging engineering problems. The selection of an appropriate FE code is dependent on the type of problem and computational cost. LS-DYNA R10.0 (Livermore Software Technology Corporation, LSTC, Livermore, CA, USA) is an explicit code developed for the dynamic analysis of non-linear problems, which requires small time steps. LS-DYNA R10.0 was found to be the most preferred choice for this kind of analysis due to the capability of solving the problems involving large deformation, easy application of SPH method and the vast variety of material models availabale for concrete and soil.
The objective of this study is to develop an efficient 3D coupled numerical model to probe the impact of pile punching on adjacent piles. The 3D coupled SPH-FEM model was generated based on the particle approximation approach and calibrated against field experiments. The established SPH-FEM model was then used to investigate the mechanism underlying pile interactions arising from the impact of pile punching.

Establishment of the 3D SPH-FEM Model for Pile Punching
SPH is a mesh-free Lagrangian method which employs a finite number of particles that carry individual mass to represent the material and form the computational domain [22]. the SPH method can be efficiently used for the simulation of dynamic problems involving large deformation due of its ability to handle large distortions by avoiding mesh tangling and remeshing, [23]. Although, SPH has great advantages in simulating many problems in engineering and science, SPH is highly expensive in terms of computation time (especially for 3D model), as a large number of small particles would be required and the time step would become very small. Thus, coupling the SPH and Lagrangian FEM mesh is a potentially good solution in overcoming element distortion, and in maintaining good computational efficiency. In this study, SPH particles are used to model the soil domain at near field, while the conventional FEM is used to model the intermediate and far-field soil medium and the piles.
In the SPH formulation, two basic steps are involved, namely kernel approximation and particle approximation. The first step is kernel approximation, where a spatial distance between particles is covered by a smooth length over which their properties are smoothed by a smoothing kernel function. The integral representation of smoothing kernel function and its derivative are described as [24], where W is the smoothing kernel function, h is the smoothing length, Ω is the problem domain and f is a field function. The commercial software LS-DYNA R10.0 was used for simulations throughout this study. It employs following cubic B-spline smoothing function, and it has been proven to be accurate and efficient [24]: Smoothing length, h, is an important parameter in the SPH method because it determines the influence area of the smoothing function, W, for each particle [24]. Since the mass of particle in SPH is assumed to be constant, the smoothing length associated with particles should vary accordingly with density. Although using variable smoothing length increase the accuracy of the results, it will increase the computational time. In this study, the smoothing length coefficient was set to be 1.05.
In the second step, being the particle approximation step, the computational domain is discretized with a set of initial distribution of particles that carry an individual mass. The field variables on a particle are estimated by a summation of the values over the nearest neighbor particles [24].
In the study, the particle approximation was used to generate the SPH-FEM model. The governing equations for SPH particles can be written as, where m is the mass, ρ is the density and v is the velocity. σ αβ is the total stress tensor, X is the spatial coordinate of the particle, t is the time, W is the smoothing kernel function and Π is the Monaghan artificial viscosity. The study first simulated and validated the experiment conducted by Nilson [25]. Nilson [25] recorded series of ground vibration measurements using the vibration sensors arranged at 10, 20 and 40 m distance from a pile drive. A reinforced concrete pile with a square cross-section of 270 mm × 270 mm and the length of 29.3 m were used in his experiment. The soil profile in the test area was 3 m of surface fill deposited on 12 m thick layer of medium stiff clay and a layer of 7 m thick sand on glacial till. Figure 1 shows the generated 3D SPH-FEM model for pile punching, which consists of piles and soils. Symmetric modelling capabilities play an important role in numerical analysis to save the computational effort [26,27]. However, in certain cases, the symmetric boundary conditions cannot be applied due to the presence of non-symmetries in loading, material and boundary conditions [28]. Considering the symmetries of the boundary conditions and applied loadings, only a quarter of the model was developed to reduce the computational cost in this study.
The domain of the soil was modelled with four different layers of soil to simulate the geotechnical soil profile at the test site and then was set to be 80 m long, 40 m wide and 40 m high. SPH particles were used to model the soils where large deformation is expected to occur near the driven pile. A preliminary analysis was carried out to determine the best size of the soil domain to model with SPH particles. It was found that numerical instabilities occurred due to a large element distortion when the domain as too small. In contrast, larger domain for SPH soil domain led to high computational cost. Higher accuracy of the analysis was ensured by using 0.5 m × 0.5 m size of the SPH soil domain around the pile. Eight-node solid elements with reduced integration and hourglass control were used to model the pile and soils in the far-field. The soil, close to the driven pile, was modeled with SPH particles and the rest of the model was modeled with the conventional Lagrangian meshes. With an equal distance of 10 mm between SPH particles at all axes, 270,000 particles were created to model the soils in the near field. The driven pile was modeled with solid elements with 25 mm edge length. The rest of the model was created using solid elements with the 250 mm mesh. The developed model has 437,090 solid elements. The nodes in the symmetry boundaries were fixed against translational displacements normal to the symmetry plane. The bottom of the mesh was modeled as fixed in all directions to prevent the boundary from moving in any direction. Non-reflecting boundaries were applied to the other surfaces, except the top surface which has the free boundary condition. A symmetry boundary was applied to those SPH particles at the symmetry planes using *BOUNDARY_SPH_SYMMETRY_PLANE.
Four different soil layers were simulated in this SPH-FEM model. Thus, the model consisted of four different SPH parts with different soil densities. Various methods exist which could handle the interactions between different SPH parts. The standard SPH interpolation functions can be used to handle the interaction between SPH parts. No contact definitions are needed and multiple SPH parts are treated as one part in the standard SPH interpolation. However, when the densities and masses of neighboring particles vary largely within the smoothing length, the standard SPH interpolation gives false values on the smoothing quantities of a particle. Muller et al. [29] showed that when the density ratio larger than 10, the interaction between SPH parts cannot be realistically simulated using the standard SPH interpolation. The instabilities, due to large density ratios across the interfaces, can be avoided by introducing a penalty based node to node contact algorithm for the interaction between two SPH parts. However, when the two SPH parts have similar density and material properties, the standard SPH interpolation method has better accuracy around the interfaces [30]. Since the soil densities do not vary significantly, the standard SPH interpolation interaction was used in the present study. To activate this, CONT parameter in *CONTROL _SPH was set to 0, and no contacts were defined between those SPH parts.
Three different methods have been involved in the coupling of SPH particles and conventional FEM meshes [26,31]. The first method is SPH particles tied to the corresponding surfaces of FEM meshes as shown in Figure 2a. If the SPH particles are not tied to the FEM mesh as shown in Figure  2b, the interaction between them is achieved by the penalty based nodes to surface contact. The third method uses hybrid elements as transit layers between SPH particles and FEM meshes as shown in Figure 2c. The tied interfaces between SPH particles and FEM elements (Figure 2a) were employed in this study to couple the soil model with SPH particles and FEM elements. The interaction between the SPH and FEM elements of the driven pile was defined using penalty based algorithms, *CONTACT_AUTOMATIC_NODES_TO_SURFACE in LS-DYNA R10.0. The slave part was defined with SPH particles and the master part was defined with finite elements (i.e., the driven pile). In this method, when a slave node is in contact with the master surface, a restoring force is applied to prevent the penetration, which is directly proportional to the penetration into the solid element. Thus, when solid elements interacted with SPH particles, the SPH-FEM coupling enabled the stress transfer at the interface without penetration of SPH particles. The restoring force, F, is defined by in Equation (6), where k is the linear spring constant, d is the penetration distance and n is the surface normal vector. The impact of the hammer was applied on the pile head as an impulse using a rectangle function for force versus time. The applied load on the top surface of the pile was derived from the mass of the hammer, m, and the height of the fall, h, as given in Equation (7), where I is the impact momentum, g is the gravitational acceleration and η is the effective ratio due to the damping of the cushion. In this study, η was taken as 0.9 for the calculations.
In this study, *MAT_CONCRETE_DAMAGE_REL3 (MAT_72R3) material model was used to model the concrete pile. The advantage of this model is that the unconfined compressive strength and density are the two parameters that are required in the automatic parameter generation to simulate the concrete behavior [4]. The concrete density, compressive strength of concrete, and Poisson's ratio were considered as 2400 kg/m 3 , 25 MPa, and 0.3, respectively. Each soil layer was modelled with *MAT_MOHR_COULOMB (MAT_173) material model and the material parameters for each soil layer are listed in Table 1.

Model Calibration
In the calibration process, the SPH-FEM model was run in two steps. The first step was stress initialization to induce steady initial in-situ gravity stresses in the soils using the *CONTROL_DYNAMIC_RELAXATION option in LS-DYNA R10.0. The impact load on the pile was then applied as the second phase after the dynamic relaxation phase. The soil-pile interactions and ground vibrations were analysed in the second phase.
Calibration of the coupled SPH-FEM modelling technique was carried out against field tests [7,25]. Massarsch and Fellenius [7] presented the results of punching one test pile obtained from a series of field test carried out by Nilson [25] in Sweden. The test pile was a reinforced concrete pile with a square cross-section of 270 mm × 270 mm. The bulk density and the impedance of the pile were 2400 kg/m 3 and 714 kNs/m, respectively. The total length of the pile was 29.3 m. The pile punching involved a 4000 kg weight hammer falling 0.4 m per blow. Ground vibrations were measured at a horizontal distance of 10, 20 and 40 m from the driven pile as shown in Figure 3. In the field test by Nilson [25], the geophones were used to measure the particle velocities vertically (V1, V2 and V3) and horizontally in the radial (H4) and transverse (H5) directions of wave propagation. Three monitoring points on the soil surface at 10, 20 and 40 m distances from the driven pile were defined using the *DATABASE_HISTORY_NODE option in LS-DYNA R10.0. LS-DYNA R10.0 offers options to extract the all nodal time history data from the nodal output. Velocity-time histories of the ground vibration at these monitoring points were extracted to compare with the experimental results. Figure 4 shows a comparison of the ground vibration results from the 3D SPH-FEM analysis and field measurements (at the pile depth of 3 m). The plots show broad agreements of the results (in terms of waveforms at the monitoring points) from the calibrated SPH-FEM model and the experiment. A common observation is that the numerical results for peak velocities at the monitoring points are slightly lower than the field test results, which is probably due to the fact that a simplified ground profile was used in the SPH-FEM model. Moreover, due to the lack of information on the hammer impact function, a rectangular function was used in the SPH-FEM analysis to apply the hammer impact on the pile head. This might be another reason for the discrepancies observed between the numerical results and field test results. Although, the results from the calibrated SPH-FEM model are somewhat lower than the field measurements, it still can be seen that the simulated results are in good agreement with the field monitoring results, which provides adequate confidence for using the established SPH-FEM model, in order to study the pile punching effects on adjacent piles.

Case Study: Impact of Pile Punching on an Adjacent Pile
The impact of pile punching on an adjacent pile was investigated using the established 3D SPH-FEM model. In the model, a 5 m long pile with 200 mm diameter circular cross-section was additionally generated ( Figure 5). The pile punching in clayey soil and sandy soil were considered in the study. The soil was idealised as homogeneous and isotropic material. Note that the influence of water table was not considered. The driven and the adjacent piles as well as soil were modelled using the same material models described in Section 2. The material properties for the clay, sand and concrete piles are given in Section 2. The interaction between the adjacent pile and surrounding soil was modelled by using AUTOMATIC_SURFACE_TO_SURFACE contact option in LS-DYNA R10.0. This assumes contact at the surface and enables transfer of stresses between solid elements. A parametric study was carried out to investigate the lateral response of an adjacent pile in clayey soil and sandy soil due to pile punching, by varying the clear spacing between the piles from 1.5 to 10d (d is the pile diameter). The impact of the hammer was applied on the pile head as an impulse using a rectangle force function with time. A load of amplitude 3285 kPa and duration 0.1 s was applied to the pile head. The period of the hammer blow was considered as 0.5 s. Figure 6 depicts the numerical results for lateral displacement of the head of the adjacent pile against the penetration depth of the driven pile when the clear spacing is 4 times the pile diameter. It can be seen from Figure 6 that the lateral displacement at the head of the adjacent pile initially increased, followed by a slight decrease as the penetration depth of the driven pile increased. Also, note that the lateral displacement of the head of the adjacent pile was smaller for a driven pile in a clayey soil compared to sandy soil. In sandy soil, the displacement of the soil was larger due to the weak bond between the soil particles. Thus, it might be the reason why the adjacent piles are expected to have a larger displacement when the piles are driven in sandy soil.  Figure 7 shows the impact of the pile spacing on the lateral displacement of the head of the adjacent pile. As expected, as the clear spacing between piles increased, lateral displacement of the adjacent pile decreased. It was also observed that the lateral pile head displacement decreased from 16 to 2 mm as the clear spacing between piles increased from 1.5 to 10d when the pile were embedded in the clayey soil. However, in sandy soil, the pile displacement decreased from 54 to 8 mm as the clear spacing between piles increased from 1.5 to 10d. Thus, the interaction between piles was considerably higher in sandy soil than that in clayey soil. A parametric sensitivity study was further carried out for different soil elastic modulus and soil density. Soil density was varied from 1600 to 2200 kg/m 3 , while the elastic modulus of soil was varied from 10 to 100 MPa to represent very soft to stiff clay soil. Figure 8 shows the lateral displacement at the head of the adjacent pile against the clear spacing between piles for different soil densities. As can be seen, the density of the soil affected the lateral displacement of the adjacent pile. As the soil density increased, the lateral displacement of the head of the adjacent pile resulted from pile punching decreased. This is because the pile is subjected to higher inertial force when it is driven in the soil which has high density. Figure 9 shows lateral displacement at the head of the adjacent pile against the clear spacing between piles for different soil elastic modulus. The elastic modulus of the soil also significantly affected the lateral displacement of the adjacent pile during the pile punching. The lateral displacement of the head of the adjacent pile caused by pile punching increased as the elastic modulus of soil decreased. Although there is some match in the initial slopes of the lateral displacement curves, it indicates a considerable difference when the pile spacing to diameter ratio greater than 3. In stiff clays, it was observed that the impact of pile punching on adjacent existing pile is comparatively less. Thus, it is clear that the interaction between piles reduces when the piles are driven in stiff soils.
Moreover, Table 2 summarises the numerical results for the lateral displacement of the head and tip of the adjacent pile when the clear spacing is 3 times the pile diameter. Tilt of the adjacent pile was calculated as the ratio of head displacement relative to the tip displacement of the pile to pile length. In all the cases, it was observed that the tilt of the adjacent pile is insignificant. The maximum tilt of 0.00304 (i.e., about 1/329) was obtained for the very loose sand which has elastic modulus of 10 MPa. The results show that tilt of pile decreases when the pile is installed in a dense or hard soil. Also, it is clear that the installation of a pile close to an existing pile will induce an overall lateral displacement of the adjacent pile rather than tilting.

Conclusions
In this study, the impact caused by pile punching on an adjacent pile was investigated using a 3D well-established SPH-FEM model; the model was calibrated against field measurements. A comprehensive parametric sensitivity study was performed to evaluate the impact of soil properties on the displacement of a pile, due to the punching of an adjacent pile, by varying the elastic modulus of the soil, soil density and spacing between piles. It was found that the lateral displacement of an adjacent pile (due to pile punching) increased with the decrease in soil elastic modulus, soil density and the spacing between the piles. The interaction between piles became weaker when the piles are driven in stiff soils. The results also show that the lateral displacement at the head of an adjacent pile was fairly significant for piles driven into sandy soil.