Implementation of spatial SV wave oblique incidence and helical pile response under this condition based on Python in finite element analysis

1.   Introduction

Recent advances have focused extensively on understanding how underground structures react when seismic waves strike at oblique angles ^[1,2]. For instance, Álamo et al. ^[3] examined how structures, piles, and surrounding soils interact, particularly considering different shear wave types and the effects of varying incidence angles. Their findings showed that oblique wave incidence has a more significant impact on these interactions compared to vertical incidence, with the wave's dominant frequency playing a crucial role. Yang and Zhou ^[4] explored how submarines respond dynamically to obliquely incident P waves and SV waves in a two-dimensional context, considering their interaction with the overlying water. Similarly, Wang et al. ^[5] introduced a technique for simulating oblique P-wave incidence in 3D using wave field superposition principles. Their analysis, carried out using the ABAQUS finite element software, delved into the dynamic responses of embankment dams with asphalt concrete cores when subjected to both horizontal and vertical oblique wave incidences. Zhao et al. ^[6] investigated the effects of obliquely incident P waves and SV waves on the interaction between structures, soil, and piles in undersea tunnels, concluding that the internal forces in tunnel linings increased with larger incidence angles. Zhou et al. ^[7] assessed the dynamic behavior of subway tunnels, the surrounding soil, and adjacent structures when impacted by oblique SV-wave incidences, factoring in aspects such as tunnel depth, spacing, and frame distance. Currently, there are few studies on fast finite element simulations of seismic wave propagation at oblique incidence in space and the effects on structures.

In the field of numerical simulation, scholars have conducted relevant studies on helical piles, but most of them are based on ordinary piles. El-Sawy ^[8] detailed the construction and validation process of the three-dimensional finite element helical pile model, and the study showed that increasing the stiffness of the pile may lead to a slight increase in the soil reaction. In-depth studies may be needed for piles of larger diameters; the addition of a second helical blade did not alter the shear or bending behavior of the original bottom spiral blade. The above findings provide an important basis for the simulation of helical blades, and the finite element model developed based on these findings can accurately simulate the behavior of helical piles under seismic conditions, which can be used as a reference for the seismic design of helical piles. Hussein and EI Naggar ^[9,10,11] investigated the response of single-vane helical piles as well as double-vane helical piles in soils through the development of an OpenSees pile-soil model and obtained the results for unsaturated soils with a single blade and a double blade helical pile. They concluded that the stress and strain of double-blade helical piles in unsaturated soil are smaller than those of single-blade helical piles. After that, they focused on the effect of different scale models on the seismic response of helical piles, established 1:1 and different scale size models through the existing experiments, took into account the saturated or not saturated soil, and then concluded that the response of soil and helical piles can be extrapolated by using the scale factor. Based on the above studies, scholars have seldom explored the response of helical piles affected by oblique incident seismic waves in a marine soft soil field.

The above-mentioned studies demonstrate that researchers have extensively analyzed the oblique incidence of seismic waves using both 2D and 3D models. However, few have focused on developing tools or plug-ins for the application of P, SV, and SH waves. Additionally, most existing research considers wave incidence in only one direction, despite the fact that seismic waves propagate in all directions in three-dimensional space. Python's scripting capabilities and high compatibility with ABAQUS give it significant advantages in this study, such as efficient automated modeling, parametric control, and batch computation capabilities. To address this, the present study applies P, SV, and SH waves in ABAQUS using Python, deriving the corresponding spatial oblique incidence equations. By simulating and constructing a helical pile structure for geotechnical sites, the influence of earthquakes on the seismic settlement of soft soil sites is studied to verify the applicability of oblique incidence theory and its implementation methods.

2.   Derivation of the theory of oblique incidence in any direction in space

2.1.   Viscos-spring artificial boundary

Viscous-spring boundaries are known for their stability and resistance to problems such as low-frequency drift and high-frequency oscillations. This method is conceptually simple, easy to apply, and widely used in structural calculations. Its general form can be expressed as ^[12]:

where $l$ refers to the artificial boundary; $t$ represents the computational time; ${\sigma _{li}}$, ${u_{li}}$, ${\dot u_{li}}$ relate to the node's stresses, displacements, and velocities along the direction of computation, while ${K_{li}}$, ${C_{li}}$ pertain to the spring and damping element parameters of the viscous-spring artificial boundary in the same direction.

Deeks and Randolph ^[13], together with Liu and Lv ^[14], introduced a two-dimensional viscoelastic artificial boundary, demonstrating that the viscous-spring boundary effectively replicates the elastic rebound in a semi-infinite medium beyond the boundary, maintaining stability across both low and high frequencies. Building on this work, Du et al. ^[15] extended the concept to three-dimensional viscous-spring boundary conditions, utilizing the empirical superposition of plane waves and scattered waves from distant fields to more accurately model the behavior of outgoing waves. The damping and spring parameters of the viscous-spring boundary can be expressed in both tangential and normal directions.

where variables $n$ and $t$ denote the normal and tangential components, while $K$ and $C$ represent the parameters for spring and damping. Furthermore, $\lambda$, $G$, and $\rho$ refer to the first Lamé constant, shear modulus, and material density, respectively. The recommended value for the stiffness adjustment factor A is 0.8, while the damping adjustment factor B is suggested to be 1.1. The propagation speeds of P waves and S waves are denoted as ${c_p} = {(\lambda + 2G/\rho )^{1/2}}$ and ${c_s} = {(G/\rho )^{1/2}}$, respectively. Additionally, the influence of viscous-spring boundary stresses on the node located at the artificial boundary is quantified by ${L_l} = \sum\nolimits_{e = 1}^N {{L_{le}}}$. Here, $N$ is the number of adjacent elements at node l, ${L_{le}}$ represents the boundary stress influence on element e, and $r$ is the distance from the node to the center of the structure.

2.2.   Equivalent seismic loading at artificial boundaries

This research utilizes an equivalent seismic loading method, employing viscous-spring artificial boundaries to convert the wave scattering problem into a wave source problem, allowing seismic waves to be introduced into the model. The general expression for this approach is given as ^[16]:

In this equation, ${F_{li}}$ denotes the equivalent nodal load, while $\sigma _{li}^F$, $u_{li}^F$, $\dot u_{li}^F$ represent the displacements, velocities, and surface stresses at the artificial boundary in response to the site conditions. Additionally, ${K_{li}}$, ${C_{li}}$ correspond to the viscous-spring and damping parameters of the viscoelastic artificial boundary, which are computed using Eqs. (2)~(5).

A method for introducing seismic waves through multiple boundary surfaces is formulated based on the approach in ^[16], enabling oblique SV-wave incidence from any spatial direction. As illustrated in Figure 1, SV waves enter the space at oblique angles from two different directions. One of the incident azimuths refers to the angle between the direction of incidence and the positive y-axis, while the other represents the angle with the positive z-axis. When the SV wave initially leaves the wavefront and reaches the two azimuths at distinct surfaces, two outcomes arise. In the first scenario, the SV wave directly impacts the surface. In the second scenario, the wave reflects off the ground, creating an SV wave with the same incidence angle and a P wave at a steeper angle upon reflection. According to Snell's law ^[17], there is a critical incidence angle for oblique SV waves. If this critical angle is surpassed, the model can become distorted, leading to inaccuracies in the calculations. This critical incidence angle is determined as follows:

Figure 1.  3D model SV-wave oblique incidence.

DownLoad: Full-Size Img PowerPoint

Initially, the delay time for the incident wave to reach the truncation plane is calculated. Specifically, $\Delta t_1^{{\text{ZP}}}$, $\Delta t_2^{{\text{ZP}}}$, $\Delta t_3^{{\text{ZP}}}$ represent the delay times for the incident SV wave, the reflected SV wave, and the reflected P wave on the ZP plane. Similarly, $\Delta t_1^{{\text{ZN}}}$, $\Delta t_2^{{\text{ZN}}}$, $\Delta t_3^{{\text{ZN}}}$ denote the delay times for the incident SV wave, reflected SV wave, and reflected P wave on the ZN plane. For the XP plane, $\Delta t_1^{{\text{XP}}}$, $\Delta t_2^{{\text{XP}}}$, $\Delta t_3^{{\text{XP}}}$ signify the delay times of the incident SV wave, reflected SV wave, and reflected P wave, while $\Delta t_1^{{\text{XN}}}$, $\Delta t_2^{{\text{XN}}}$, $\Delta t_3^{{\text{XN}}}$ correspond to these same delay times on the XN plane. Lastly, $\Delta t_1^{{\text{YP}}}$ indicates the delay time for the incident SV wave on the YP plane, as shown below:

where ${L_x}$, ${L_y}$, ${L_z}$ are the lengths of the truncated body on the x, y, z axes, respectively, and x, y, z are the coordinates of point A. The displacement of the ZP surface consists of the incident wave SV wave displacement ${u_1}(t - \Delta t_1^{{\text{ZP}}})$, the reflected wave SV wave ${u_2}(t - \Delta t_2^{{\text{ZP}}})$, and the reflected wave P wave ${u_3}(t - \Delta t_3^{{\text{ZP}}})$, as follows:

where $u_x^{{\text{ZP}}}$, $u_y^{{\text{ZP}}}$, $u_z^{{\text{ZP}}}$ are the displacements in the x, y, and z directions; ${B_1}$ and ${B_2}$ are the amplitude ratios of the reflected SV wave and the reflected P wave to the SV wave, respectively, as shown below ^[18]:

The displacement equations on the three surfaces ZN, XP, and XN are, respectively:

The YP surface displacement equation is:

The stresses at each equivalent node can be calculated at the ZP surface based on the free-field displacement equation:

where ${\dot u_1}(t - \Delta t_1^{{\text{ZP}}})$, ${\dot u_2}(t - \Delta t_2^{{\text{ZP}}})$, ${\dot u_3}(t - \Delta t_3^{{\text{ZP}}})$ are the displacements in the x, y, and z directions, respectively. Similarly, the stress equations on the ZN, XP, XN, and YP surfaces can be calculated:

Finally, only the displacement, velocity, and stress equations on the corresponding surfaces need to be substituted into Eq. (6) to obtain the corresponding equivalent nodal loads.

2.3.   SV wave input validation

In this study, the viscous-spring artificial boundary and the corresponding equivalent nodal forces were implemented in ABAQUS using a custom Python script, as shown in Figure 2. To assess the accuracy of the equivalent nodal load, this section simulates the propagation of SV waves in arbitrary space through numerical modeling. The truncated body's dimensions are set to 1200 m × 600 m × 200 m, with an elastic modulus of 5 GPa, a Poisson's ratio of 0.3, and a mass density of 2500 kg/m³. The observation point is located at point A (x, y, z = 600 m, 600 m, 200 m). A pulse wave, shown in Figure 3, is chosen as the input, and its displacement-time curve is defined as follows:

where $G = {\left( {\frac{t}{{{T_0}}}} \right)^3} \cdot H\left( {\frac{t}{{{T_0}}}} \right)$, H is the Heaviside function; the peak displacement ${u_0}$ is taken to be 1 m, and the duration of the pulse wave action ${T_0}$ is taken to be 0.25 s.

Figure 2.  Flowchart of SV-wave imposition in ABAQUS.

DownLoad: Full-Size Img PowerPoint

Figure 3.  Displacement history of the incident SV wave.

DownLoad: Full-Size Img PowerPoint

2.4.   Boundary conditions and analysis steps

Since the conventional method of applying viscoelastic artificial boundaries and equivalent seismic loads has the problems of a slow application process and being error-prone, a plug-in has been written through Python that can quickly add viscoelastic artificial boundaries and equivalent seismic loads, and secondly this plug-in can also perform geostress balancing and dynamic stress balancing. For the research needs of the seismic dynamic response of the structure under the influence of oblique incidence, based on ABAQUS finite element software and the Python program, combined with the viscoelastic boundary and equivalent seismic load input technology, a plug-in program is developed, which contains the module of geostress equilibrium and the module of seismic analysis. The module can automatically complete the equilibrium of geostress under static boundary conditions and dynamic boundary conditions (dynamic and static boundary conversion equilibrium must be completed, otherwise the model does not conform to the actual situation and cannot be calculated); the seismic analysis module can output the information of the boundary nodes of the model, apply the viscoelastic boundaries, and automatically calculate and add the equivalent seismic load at each boundary node location. In this part of the study, the plug-in shown in Figure 4 is used to apply the viscous-spring artificial boundary and the equivalent seismic loads. The initial stresses, displacements in the structure, and the original mass of the soil under the influence of gravity are handled as follows:

Figure 4.  Python plug-in for finite element simulation software.

DownLoad: Full-Size Img PowerPoint

(1) Step 1: A static analysis is carried out, taking into consideration the gravitational effects on the soil. The earth stress equilibrium is achieved using the ODB (object data base) method, and the process is completed when the soil displacement reaches 10⁻⁴ m. Since an oscillatory response may arise during the transition from static to dynamic boundaries ^[19], this issue is addressed, and the initial conditions are determined using dynamic techniques.

(2) Step 2: This step involves dynamic analysis, which is executed once the equivalent seismic load has been applied.

The plug-in program has a high degree of integration, fully interfaced operation, no need for programming, and is convenient for application in practical engineering. It can automatically complete the tedious work such as geostress equilibrium iteration and the calculation of viscoelastic boundary parameters and equivalent seismic loads, which greatly improves the efficiency of model pre-processing work.

Figure 5 illustrates the displacement cloud at various time intervals when the SV wave obliquely strikes the truncated body at angles of $\alpha$ = 10, 20, and 25°, $\beta$ = 10, 20, and 25°. It is observed that when the incident wave does not interact with the ground, no reflected waves are generated. However, once the wave reaches the ground, it splits into an incident SV wave, a reflected SV wave, and a reflected P wave due to ground reflection. As shown in Figure 5, the incident wave forms angles of 15°, 20°, and 25° with the xz-plane, demonstrating the success of the SV-wave oblique incidence simulation in the horizontal plane.

Figure 5.  Displacement propagation of space body SV waves at different times for $\alpha$ = 15, 20, and 25°, $\beta$ = 10, 20, and 25°.

DownLoad: Full-Size Img PowerPoint

Figure 6 displays the time-range curves for horizontal and vertical displacements at point A, located at the center of the surface, for incidence angles of $\alpha$ = 10, 20, and 25°, $\beta$ = 10, 20, and 25°, compared with the analytical solution. The numerical results align well with the analytical solution, validating the accuracy of the spatial SV-wave oblique incidence input direction proposed in this paper. In addition, its amplitude at the surface is greater compared to the incident wave amplitude, which is mainly due to the fact that the seismic wave is incident from the bottom of the soil layer, and when it reaches the top of the soil layer, its acceleration undergoes a method, and is therefore not the same as at the time of the incident.

Figure 6.  Comparison of horizontal and vertical displacement curves with the analytical solution at point A.

DownLoad: Full-Size Img PowerPoint

3.   Simulation case of helical pile seismic settlement in the geotechnical field

3.1.   Finite element model

To study the vertical displacement, seismic subsidence, and the dynamic p-y curve of the helical pile during the oblique incidence of the P wave, a three-dimensional finite element model of the helical pile-soil system is initially established. This model primarily consists of the soil and the helical pile which is located in the central region of the soil and has an embedment depth of 3 m, as shown in Figure 7. The direction of the x-axis is the incident direction of the seismic wave, and the y-axis is the vertical direction. Since there is an inclination angle in the helical blade of the helical pile in the actual project, the helical blade is made into a plane disc to simplify the model, and this reduces the calculation error compared to the case of having an inclination angle, and the influence of the inclination angle can be neglected ^[20].

Figure 7.  Helical pile-soil finite element model.

DownLoad: Full-Size Img PowerPoint

For seismic waves to propagate efficiently through the soil, soil elements must conform to a specific range of dimensions to ensure the computational efficiency and reliability of the model ^[21]:

where ${V_s}$ is the shear wave velocity of the soil layer and ${f_{\max }}$ is the cutoff frequency, which is generally selected as the highest frequency of the ground-shaking band. Given that the input seismic wave predominant frequency in the model is 2.5 Hz and the minimum shear wave velocity is 71 m/s, substituting this frequency into Eq. (25) results in ${h_{\max }} \leqslant$ 4.73 m. The maximum element size in the model is 0.15 m, satisfying the accuracy requirements. Finally, the soil is divided into 29,196 elements.

To reduce the influence of soil property changes in the buried area on the seismic subsidence of the helical pile, it is necessary to construct the seismic subsidence analysis model of the homogenized soil layer. Moreover, the interpolation method is used to homogenize the soil layer within a 5-m depth, as shown in Table 1. The soil uses an elastoplastic constitutive model with Mohr-Coulomb as the damage criterion, which responds excellently to the mechanical properties of soft soil. The two main research objects in this paper are the helical pile seismic subsidence and the dynamic p-y curve. In the case of the dynamic p-y curve, p is the reaction force of the soil on the side of the pile on the pile, which is a lateral load, and y is the relative displacement of the pile and the soil in the far field on the side of the pile ^[22]. Previous studies have indicated that the load-displacement response of transversely loaded piles is unaffected by the Mohr-Coulomb model ^[23]. Consequently, numerous researchers have employed the Mohr-Coulomb model as a means to investigate the damage criterion for laterally loaded piles ^[24,25]. In the case of seismic subsidence, scholars have developed elastic-plastic constitutive models incorporating the Mohr-Coulomb damage criterion, utilizing the softened model analysis method ^{[26,27,28,29,30]}.

Table 1.  Soil constitutive model and contact model parameter table.

Material	Density/kg·m-3	Young's modulus/MPa	Poisson's ratio	Friction angle/°	Cohesive force/kPa	Coefficient of friction
Soil	1890	2	0.30	15	15	0.1

---

 | Show Table

DownLoad: CSV

3.2.   Seismic wave

Two distinct seismic waves, namely the El-Centro and Kobe waves, have been selected for this study to simulate the seismic response of helical piles across a spectrum of seismic frequencies. The peak accelerations of the two seismic waves were adjusted to 0.05 g, 0.1 g, and 0.3 g, respectively, to accurately simulate the seismic subsidence of the helical piles under varying PGAs (peak ground acceleration). The acceleration time history curves and the corresponding Fourier spectra of the two seismic waves are presented in Figure 8. The frequency range of the El-Centro wave is approximately 0.5 Hz, while that of the Kobe wave is approximately 1 Hz. To investigate the impact of various incidence angles on the seismic response of helical piles, six distinct angles were considered: 0°, 15°, 30°, 45°, 60°, and 75°. As outlined in the previous section, the influence of seismic subsidence on helical piles was investigated, considering various factors such as incidence angles, seismic waves, and PGAs in soft ground.

Figure 8.  El-Centro and Kobe wave acceleration time history curves and Fourier spectral curves.

DownLoad: Full-Size Img PowerPoint

3.3.   Model validation

In order to verify the applicability of this model, the same pile-soil model is established according to the case given in the literature ^[10]. The maximum horizontal acceleration and soil displacement are analyzed, shown in Figure 9. It can be seen that the model established in this section can fit the acceleration and soil displacement curves given in the literature better and the error value is basically within 10%. Although the acceleration error at the pile head has exceeded 10%, this value is on the high side, which adds safety and reliability to the seismic design of the pile. Therefore, the method proposed in this paper is suitable for analyzing the wave input problem of the pile-soil model.

Figure 9.  Validation of the pile-soil model.

DownLoad: Full-Size Img PowerPoint

3.4.   Response results of the helical pile model

3.4.1.   Effects of different seismic waves

Figure 10 depicts the vertical displacement and horizontal displacement within the xy cross-section of the double-helical blade pile-soil model subjected to various seismic waves. The El-Centro wave is observed to induce greater vertical displacement in the pile-soil model than the Kobe wave. Owing to the differing frequencies of the seismic waves, the helical pile-soil model demonstrates heightened sensitivity to the El-Centro wave. Furthermore, the displacement cloud of the model shows a predominance of the blue spectrum under El-Centro waves, whereas a red spectrum is more prevalent when subjected to Kobe waves. This indicates that under El-Centro waves, numerous regions exhibit larger vertical displacements than those under Kobe waves. Therefore, considering the impact of seismic wave frequencies on the vertical displacement of the helical pile-soil model is essential. It can be seen by the horizontal displacement curve on the right side that the upper part of the deformation of the helical pile is the largest, which reaches about 0.8 mm, and the horizontal displacement is very small in the burial depth below 1 m. The horizontal displacement of the helical pile is about 0.8 mm. It can be found that the incidence angle degree has little effect on the horizontal deformation of the pile, so only the vertical deformation of the helical pile is studied subsequently.

Figure 10.  Effect of different seismic waves on vertical displacement and horizontal displacement of helical pile-soil.

DownLoad: Full-Size Img PowerPoint

3.4.2.   Effects of different peak accelerations

Figure 11 presents the vertical displacement cloud of the xy cross-section for the two-blade helical pile-soil model under various PGA conditions. The maximum vertical displacement across all pile-soil models occurs at a PGA of 0.3 g, while the minimum is observed at 0.05 g. This indicates a direct relationship between PGA and the vertical displacement of the helical pile-soil model. Furthermore, with increasing PGA, both the spatial extent of soil effects surrounding the helical pile and the magnitude of vertical displacement correspondingly increase. In this context, the helical pile's influence range in the blade region is significantly larger than in regions without blades. This phenomenon is attributed to the helical pile's downward displacement, creating a specific soil contact area in regions lacking blades. Consequently, this limits the influence range, resulting in localized soil effects. Conversely, the presence of blades expands the soil contact area, thus enlarging the influence range. This, in turn, reduces the vertical displacement experienced by the helical pile.

Figure 11.  Effect of different seismic PGAs on the vertical displacement of helical pile-soil.

DownLoad: Full-Size Img PowerPoint

3.4.3.   Effect of different angles of incidence

Figure 12 depicts the vertical displacement response of the two-blade pile-soil model at various incident angles. It is observed that the vertical displacement of the helical pile-soil model varies with changes in the incident angle. At an incident angle of 0°, the seismic wave's effect on vertical displacement is minimal. Conversely, at an incident angle of 15°, the vertical displacement reaches its peak among the six angles studied. Subsequently, as the incident angle increases, vertical displacement decreases. This demonstrates that the incident angle significantly impacts pile-soil vertical displacement, underscoring the importance of considering incident angle effects in such studies.

Figure 12.  Effect of different incident angles on the vertical displacement of helical piles.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

In this study, the viscous-spring artificial boundary is utilized, and obliquely incident angles from two different directions are considered. The corresponding equivalent nodal force equations for the spatial oblique incidence of SV waves are derived. Building on this, the oblique incidence of various types of seismic waves in a soil body is simulated using a Python program. The key conclusions are as follows:

(1) Equivalent nodal force equations for simulating oblique incidence in 3D space using finite element methods were derived, offering a more precise approach for seismic wave loading.

(2) The Python-based programs for the viscous-spring artificial boundary and equivalent seismic load inputs addressed challenges such as complexity in ABAQUS, high computational workload, and difficulty in identifying errors, thereby ensuring the reliability and practicality of the numerical solution.

(3) The accuracy of the theoretical derivations and finite element implementation for 3D spatial oblique incidence of SV waves was confirmed through comparison with analytical results.

(4) By simulating the seismic response of helical piles under oblique incidence in soft soil fields, it is demonstrated that the influence of oblique incidence on the internal structure of the geotechnical field cannot be ignored, further verifying the applicability of the method proposed in this paper.

Author contributions

Conceptualization, D.H.; Methodology, H.W.; Software, Z.H.; Validation, D.H.; Formal analysis, C.X., and S.Y.; Investigation, H.C.; Resources, D.H.; Data curation, W.Z.; Writing—original draft preparation, H.C.; Writing—review and editing, D.H., and H.C.; Visualization, H.C.; Supervision, Z.Z.; Project administration, D.H.; Funding acquisition, D.H. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (No. 52408537), the Jiangsu Province Key R & D Program (Social Development) Project of China (No. BE2021681), the Natural Science Foundation of Jiangsu Province (No. BK20230693), the China Postdoctoral Science Foundation (No. 2024M752370), the Natural Science Foundation of the Jiangsu Higher Education Institutions of China (No. 23KJB560005), the Key R & D (Social Development) Project of Lianyungang City (No. SF2327), Lianyungang City Double Creation Doctoral Program, and Jiangsu Province Graduate Practice Innovation Program Project (No. SJCX24_2078).

Conflict of interest

All authors declare no conflicts of interest in this paper.