Simplified computational methods for estimating dynamic impedance of batter pile groups in homogeneous soil

A simplification of existing analytical solutions together with an extended hybrid model is proposed for estimating the horizontal impedance of batter pile groups subjected to harmonic loading. The solutions are given through compact interaction factors and analogue pile‐soil‐pile interaction elements. Benefits and pitfalls of the simplified models are elucidated through a comprehensive comparison with a refined finite element model. The effect of varying center‐to‐center pile distance with respect to depth is discussed. It is shown that the analytical model is able to represent pile‐soil‐pile interaction rather well within the specified framework.

Easy access to commercial FE-software allows for accurate assessment of most structural and geotechnical problems. However, numerical solutions of dynamic soil-pile-structure response are both time-consuming and complex in nature, often demanding cross-disciplinary set of skills. Therefore, practical engineering often favours simple analytical or analogue methods that are robust, easy to use and verifiable. Due to the above-mentioned discouragement of inclined piles in seismic areas, researchers have in the last decades mainly focused on vertical pile groups in the development of computational methods. Kaynia [14] proposed a meticulous solution elucidating dynamic pile-soil-pile interaction and validating the superposition principle for dynamic response. Dobry and Gazetas [15] proposed a simple method for estimating the dynamic impedance of a pile group by directly applying simple wave attenuation functions as interaction factors. Gazetas and Makris [16,17] further studied the interaction factors in a two-part article series, where the inertial and flexural resistance of the receiver pile was recognized for lateral interaction. Makris and Gazetas [18] also investigated the effect of phase differences in the interaction factor approach. Mylonakis and Gazetas [19] extended the interaction method by taking into account finite pile length and soil layering. Takewaki and [20] applied the principles of the interaction factor method in order to estimate the interstory drifts in buildings. Wang et al. [21] extended the interaction factor method by including shear deformations and rotational inertia of the piles and shear deformations of soil. In the recent years however, attempts have been made to develop simplified methods also for batter pile groups. Ghasemzadeh and Alibeikloo [22] presented a simple, closed-form solution for infinitely long piles in homogeneous soil. A similar approach was presented by Ghazavi et al. [23] . Wang et al. [24] extended his shear and multi-layer model to inclined piles. Goit and Saitoh [25] performed an experimental study to asses the impact of non-linearity on interaction factors for batter pile groups and proposed additional multiplication factors.
The objective of this paper is (1) to present simple, spreadsheet friendly expressions for the dynamic impedance of batter pile groups in linear, homogeneous soil, (2) contribute in the development and understanding of simplified computational methods for pile groups and (3) illuminate the potential benefits and pitfalls intrinsic to analogue simplifications. First, we will present a simplified, closed-form solution that closely follows the methods for dynamic pile-soil-pile interaction available in the literature. [15][16][17][18]22,23] The presented method is valid for long, floating piles with fixed head conditions in homogeneous soil. The underlying philosophy in approaching this method is primarily motivated by (1) simplicity in use, (2) accuracy in conjunction with the intended field of application and (3) educational value. Next, the analytical solution is evaluated by computing the lateral, dynamic impedance of a two-by-one pile group. The effects of pile distance, batter angle and load frequency are evaluated and compared to a fully numerical model using beams, solids and appropriate boundary elements constructed in OpenSees MP. [26] Finally, a hybrid model is constructed where the piles are modelled as finite element beams, while the soil and the pile-soil-pile interaction is represented using analogue elements based on the analytical approach. The hybrid model may be considered as an extension of the analytical model, where the limitations of the closed-form solution are, at least partly, bested by solving the differential equation numerically.

Simplified analytical model
The general aim is to establish expressions for the displacement of a receiver pile due to harmonic vibrations of a neighboring pile through an interaction factor. The main principles for this approach are well-established throughout the literature [14,15,[19][20][21][22]24] and will not be repeated here. One of the more crucial assumptions is that wave emission occurs simultaneously from all points along the pile. This assumption is thoroughly discussed by Makris and Gazetas [18] who have shown that the assumption is reasonable for vertical piles. For batter piles, even though the wave emission from the source pile may be regarded as simultaneous, the waves striking the receiver pile are generally not in phase along the pile length due to the batter angle. This issue is discussed in Section 2.1.4. Past studies have shown that rotational interaction is negligible for vertical piles. [14,27] The same arguments apply to batter piles. The limitations of the method presented in this paper include inertial loading, long (floating) piles, homogeneous deposits, cylindrical pile shapes and fixed-head conditions. Note that other boundary conditions may of course be enforced, but every specific combination of conditions requires an unique solution. The piles are modelled as bars (axial loads) and Euler-Bernoulli beams (lateral loads). The soil is modelled using linear, frequency-dependent Winkler springs and dashpots. [18,28,29] A schematic illustration of the model is shown in Figure 1. Expressions used for the springs are ≃ 0.6 F I G U R E 1 Schematic illustration of a harmonically loaded two-pile group where subscripts and denote axial and lateral deformation, respectively. is the Young's modulus of soil and 0 is the dimensionless frequency where is pile diameter, is the angular excitation frequency and is the shear wave velocity. The expressions used for dashpots consist of radiation damping and hysteric damping, that is, Here, is the soil density, is hysteric damping factor of the soil (not to be confused with the batter angle in Figure 1) and is the Lysmer's analogue velocity. The displacement vector for a pile group may be expressed as where (axial) and (lateral) are the total receiver pile displacements, (axial) and (lateral) are the receiver pile displacements due to own loading, (axial) and (lateral) are the source pile displacement due to own loading, is the axial-axial interaction factor, is the axial-lateral interaction factor, is the lateral-axial interaction factor and is the lateral-lateral interaction factor. Due to the rigidity of the pile cap, the total displacement of the single pile is approximately equal to the pile group displacement.

2.1.1
Axial response of a single pile Assuming harmonic, steady-state response, the equilibrium equation of the axially vibrating pile is where is the E-modulus of the pile and is the cross sectional area of the pile. The general solution of the free vibration response is The roots are where Applying the conditions (0, ) = (0, ) = 0 and insuring a finite displacement as z tends to infinity, 1 must be equal to zero. The displacement is then expressed as 2.1.2 Lateral response of a single pile The governing differential equation is expressed as where is the second moment of area of the pile cross section. Solutions are sought separately for < and > , where For simplicity, derivations are only presented for < since this solution covers the frequency range of practical interest. The solution for > is obtained through a nearly identical procedure. The solution of Equation 10 is obtained by applying the Laplace transformation while directly incorporating the boundary condition (0) = 0 (no rotation at the pile head). To the best of the authors knowledge, this particular technique was first applied by Makris and Gazetas [17] in relation to solving the beam-on-complex-springs differential equation. The procedure yields, algebraically at least, a fairly simple solution to Equation 10, that is, where 0 is the lateral pile head displacement and and = cos 4 + sin 4 > 0, = cos 4 − sin 4 > 0 The full derivation is somewhat tedious and is omitted due to brevity.

Attenuation functions
The use of simplified plane-strain attenuation functions is well-established and discussed throughout the literature [15][16][17][18][22][23][24] and is therefore briefly addressed in this paper. When the pile is vibrating in the axial direction, cylindrical SV-waves emanate from the pile surface. The attenuation function is then given as It is assumed that the laterally vibrating pile emanates P-waves in the direction of displacement and SH-waves perpendicular do the direction of displacement. [15][16][17][18] The attenuation functions are thus where we have introduced as the angle between the direction of vibration and the center-to-center line between the two piles. The wave field for any value of between 0 • and 90 • may be expressed as with sufficient accuracy. [15,17] However, it should be noted that when the pile spacings are very short, the scattering effect will inherently reduce the interaction. Therefore, the proposed method is applicable for spacings approximately equal to or greater than 3 .

2.1.4
Axial displacement of receiver pile due to axial displacement of source pile - The flexibility of the receiver pile will resist the soil deflection, resulting in a modified displacement profile. The batter angles are, by definition, positive clock-wise. If damping is neglected when deriving the axial response of a single pile, Equation 9 is simplified to Figure 2a shows the normalized axial displacement plotted against the ratio between depth and pile diameter at =0 using Equations 9 and 19. For the frequency range of interest, the results show that even though the actual wave propagating down the pile is completely different for the damped and undamped system, the soil displacement distribution is very similar. For higher frequency values however, the displacement distribution is somewhat different. In order to simplify the algebraic expressions, Equation 19 will be used in the following derivations and examples. A similar simplification was proposed by Makris and Gazetas [17] for lateral response. The mathematical formulation for the displacement of the receiver pile is the non-homogeneous differential equation for the dynamic bar. The bar is loaded with a distributed load F I G U R E 2 Normalized displacement of source pile with and without damping plotted against the ratio between depth and pile diameter at time =0. ∕ =40, / = 1000, damping ratio = 0.05, = 0.75 , and = 0.25

F I G U R E 3 Attenuation function amplitude.
= 100 m/s, damping ratio = 0.05, = 0.75 and = 0.25 equal to the induced soil displacement multiplied by the soil impedance, that is, For vertical piles, the center-to-center distance ( ) between two piles is constant along the pile length. For inclined piles however, the distance S varies linearly. With reference to Figure 1, the distance between two piles at an arbitrary depth measured radially from the vibrating pile is expressed as By introducing the cardinal variable in the attenuation function, the solution of Equation 21 becomes inconveniently complex. In order to investigate the effect of a varying distance ( ), Equations 15 and 16 are further explored. The first and second factor express the decay in displacement amplitude due to radiation and hysteric damping, respectively. The third term expresses the phase lag at a distance ( ). Figures 3 and 4 show the total amplitude and phase lag as functions of ( ). Firstly, it is observed that the attenuation function amplitude decays rapidly in vicinity of the source, but that the function smoothens out as the distance increases. The observation is independent of frequency. Secondly, the phase lag becomes quite prominent for both shear and pressure waves as the frequency increases. When the piles are vertical, ( ) is constant and the phase lag is solely governed by the finite wave velocity of waves propagating down the source pile. As previously mentioned, this assumption is thoroughly discussed by Makris and Gazetas. [18] For batter piles however, the phase lag contains an additional factor introduced by the varying attenuation function along the pile length. In addition, the interaction in the vicinity of the pile cap is further complicated and the plane-strain attenuation functions together with the mathematical and geometrical treatment of the loaded receiver pile are likely to be less valid. As stated before, applying a varying distance ( ) yields an exceedingly complex solution. Since we are seeking to develop simple approximate models, such solutions would arguably contribute to the opposite, especially if applying a constant attenuation function results in errors of acceptable magnitude. We will therefore proceed the derivations by assuming a constant attention function, and assess the possible shortcomings in a finite-element and hybrid model comparison. Assuming a constant distance ( ) = (0), the general solution is obtained as In order to determine the constants in Equation 23, boundary conditions are enforced. First, it is observed that must vanish if ( ) is to remain finite as approaches infinity. Second, the normal force at the top of the receiver pile must be zero. The two boundary conditions give the total solution The interaction factor is established by dividing Equation 25 by Equation 19. At the pile head ( =0), the interaction factor for inertial loading is F I G U R E 5 Dynamic axial-axial interaction factor for batter piles. ∕ =40, / = 1000, damping ratio = 0.05, = 0.75 , and = 0.25 Figure 5 shows that the axial-axial interaction decreases as the batter angles increases. The use of a constant attenuation function for axial-axial interaction may thus be argued on the basis of (1) the decreasing displacements along the pile length and (2) the fact that axial-axial interaction decreases with increasing batter angle.

2.1.5
Lateral displacement of receiver pile due to axial displacement of source pile - The soil displacement generated from the axially loaded batter pile also induces lateral displacements at the receiver pile. The flexibility of the beam will resist the soil deflection, resulting in a modified displacement profile. The mathematical formulation for the displacement of the receiver pile is the non-homogeneous differential equation for the dynamic beam. The beam is now loaded with a distributed load equal to the induced soil displacement multiplied by the soil impedance, that is, If we assume a constant attenuation function, the general solution is obtained as and In order to determine the constants in Equation 28, boundary conditions are enforced. First, it is observed that and must vanish if ( ) is to approach zero as approaches infinity. Second, the slope of deflection and the shear force are both zero at the pile head. The two boundary conditions yield  Figure 6A compares the total, particular and homogeneous solution presented in Equation 28 with constants determined from Equation 31. It is clear that the homogeneous part is negligible except for a small contribution at the top. Therefore, the homogeneous part may for all practical purposes be neglected. The lateral pile deflection may thus be reduced to The interaction factor is obtained by dividing Equation 32 by Equation 19. At the pile head (z=0), the interaction factor is Figure 7 shows that the lateral-axial interaction increases as the batter angles increase, which argues against the use of a constant attenuation function.  Figure 2B shows the normalized displacement plotted against the ratio between depth and pile diameter at = 0 using Equations 12 and 34. It is clear that the soil displacement distribution is very similar within the frequency range of interest.
In order to simply the algebraic expressions, Equation 34 will be used in the following derivations and examples. The mathematical formulation for the lateral displacement of the receiver pile due to lateral displacement of the source pile is the non-homogeneous differential equation for the dynamic beam, that is, Regarding the varying attenuation function depicted in Figures 3 and 4, the same observation are made for lateral vibrations as for axial vibrations. We will proceed with a constant attenuation function due to the substantial simplification of the algebraic expressions obtained by solving the governing differential equation. The Laplace method is yet again applied in order to obtain algebraically simple expressions. Using the same procedure as for the single pile response, the total solution may expressed as The interaction factor is expressed by dividing Equation 36 by Equation 34. At the pile head ( =0), the interaction factor is Note that when 1 = 2 = 0, Equation 37 reduces to the interaction factor presented by Makris and Gazetas [17] for vertical piles. Figure 8 shows that the lateral-lateral interaction decreases as the batter angles increases. The use of a constant attenuation function for lateral-lateral interaction is arguably reasonable do to the decreasing displacements along the pile length and the fact that lateral-lateral interaction decreases with increasing batter angle.

2.1.7
Axial displacement of receiver pile due to lateral displacement of source pile -Finally, the mathematical formulation for the axial displacement of the receiver pile due to lateral displacement of the source pile is Assuming a constant attenuation function, the general solution is obtained as The boundary conditions imply that must vanish if ( ) is to remain finite as approaches infinity. Also, the axial force at the top of the receiver pile must be zero. Enforcing the boundary conditions, the total solution may be expressed as . It is clear that the particular solution is negligible except for a small contribution at the top. Therefore, the particular solution may for all practical purposes be neglected. Note that the differential equation is a function of the spatial coordinate . Ergo, the homogeneous part is present throughout the harmonic, steady-state motion. The axial receiver pile deflection may thus be reduced to (43) Figure 9 shows that the axial-lateral interaction increases as the batter angles increases, which argues against the use of a constant attenuation function. However, the axial-lateral interaction factor is generally small compared to the other interaction factors.

Finite element model
The finite element model is constructed in OpenSees MP [26] together with the pre-and postprocessing tool STKO. [30] The piles are modelled with beam elements and the soil is modelled with eight-noded hexahedral elements. Elastic constitutive models are applied for both beams and solids. The interaction between beams and solids is modelled through rigid-linkconstraints connecting each beam node to the corresponding soil nodes such that the pile section in the given beam node acts like a rigid disk. The constraints are enforced using penalty functions. [31,32] The penalty values are obtained by approximating the order of largest entry of the stiffness matrix, that is, See Figure 10 for an illustration. The entire model is a half-space, that is, it is symmetric about the global y-axis. The geometrical cross section parameters of the beam are thus applied manually based on the half-circle and half the intended load is applied. The half-space measures 82 × 21 × 60 . The outer boundaries are represented using viscous zerolength elements developed by Lysmerand and Kuhlemeyer [33] where the damping coefficients are determined as the product of soil density, wave velocity and boundary area. Pressure wave velocity is used for damping coefficients normal to the boundary, while shear wave velocity is used for damping coefficients tangential to the boundary, that is, Here, is the boundary area represented by the given node. The zero-length elements are connected to fixed nodes. Spurious damping is avoided by choosing sufficiently smooth element mesh. The maximum element size is determined based on the smallest wave length, that is, where and 0, are the highest frequency values applied in the analysis. In order to correctly represent a wave, it necessary to use approximately ten nodes per wavelength. [34] Thus, the element size is kept below approximately ℎ = , 9 (48) for elements that seemingly contribute to the interaction. It should be noted that the element size is substantially smaller in the vicinity of the pile cap where most of the interaction is expected to occur. The final element mesh is established by verifying that refinement does not alter the results significantly. Rayleigh-damping is applied by choosing the load frequency as the first frequency. Each beam is subjected to a lateral point load applied at the top of the beam. The system is solved using Newmark's method with =

Description of analysis and model verification
The analytical method is evaluated and compared to the finite element model by calculating the dynamic, lateral impedance of a two-by-one pile group. The results are given as normalized values by dividing the actual stiffness with the summed stiffness of two vertical, individual piles. Pile spacing of 3, 5 and 10 times the pile diameter are considered together with batter angles of 0, 5, 10, and 15 degrees and homogeneous stratum with ∕ -ratios of 200 and 1000. There has been performed a total of 288 three-dimensional finite element analyzes including single pile analyzes for each batter angle and ∕ -ratio. The chosen batter angles are considered within the realistic range in practical design. The piles are free-floating with fixed-head conditions and have a length to diameter ratio ∕ = 40. A horizontal, harmonic load is applied at top of each beam. The analytical solution is derived based on the assumption that free-vibration response is negligible. Except for higher load frequencies in some system configurations, this assumption is indeed confirmed by both the finite element and hybrid model. Even when there is a notable contribution, it can be shown that the free-vibration response quickly dies out. Also, the evaluation of dynamic impedance, rather than for the maximum displacement value, pertains intrinsically to the steady-state vibration. For the hybrid and finite element model, the load is applied for a duration of 1.00 seconds using 200 time steps, and the stiffness values are obtained by averaging the absolute value of maximum and minimum displacements during the steady-state response. The models are verified by comparison to the rigorous model developed by Kaynia. [14] Figure 11 shows the real and imaginary parts and Figure 12 shows the absolute values. The models are in acceptable agreement.

Results
The results are shown in Figures 19 -18. Overall, the analytical model yields satisfactory results. The trends affiliated with pile distance and frequency are captured particularly well. In some cases, the analytical model for vertical piles tends to overestimate the group stiffness. This can be seen in Figures 14A, 16A, 17A, and 18A. The overestimation is more prominent for the higher / -ratio. The results for the 15-degree configuration with / = 200 presented in Figures 19D, 14A, and 15D shows agreement between the two models. For the 15-degree configuration with / = 1000; however, the analytical method tends to underestimate the group stiffness. The deviation increases with increasing pile distance and frequency. The mismatch is particularly notable in Figure 18D. The above-mentioned observations may be viewed from a different perceptive. Studying the results for each pile distance and ∕ -ratio, it is seen that the analytical model is somewhat less sensitive to batter angle compared to the finite element model. Generally, an overestimation for vertical piles shifts to an underestimation as the batter angle increases. This observation might indicate limitations regarding the simplified treatment of (1) wave emission and (2) how the receiver pile is loaded when the pile axis are not aligned.
Taking into consideration that the complexity of the problem is primarily governed by the pile-soil-pile interaction, it is crucial to investigate the accuracy of the analytical model in context with the dynamic response of single piles. It is observed that the discrepancy between the analytical and the finite element model in some cases roughly equals the discrepancy between the pile-group and single pile response. For example, Figure 18d shows that the analytical model better matches the single pile finite element model than the two-by-one pile finite element model for higher frequencies.
Then, the obvious argument is that if we in fact are seeking approximate results, perhaps interaction may be neglected all together. Even if that may be argued for some cases, the results clearly show that for a majority of configurations, the analytical model does in fact capture the dynamic pile group effect. It should also be added that the results in this study show that the two-by-one pile group is substantially more sensitive to load frequency compared to the single pile. The model yields reasonable results also for larger pile groups. Figure 20 shows a comparison between the proposed model and a FEM model for a three-by-three pile group consisting of both vertical and batter (15-degree) piles. Similar to the two-by-one pile group, the results show that the trends are generality well captured with reasonable accuracy.
The results presented in this study further imply that the pile-soil-pile interaction is a complex engineering problem which is not easily simplified for the general state. Indeed, the analytical model is able produce the trends affiliated with pile distance and frequency rather well, but the general inaccuracy of such methods must be recognized. In addition to the inaccuracy introduced by the various simplifications, linear methods are strongly limited in the sense that they cannot capture material (soil) and geometrical (formation of gap) non-linearity, both of which are to be expected during strong earthquake shaking and corresponding inertial loading.
The method is also limited to uniform soil profiles. Soil non-uniformity in combination with batter piles yields complex wave propagation in addition to cumbersome solutions of both the homogeneous and non-homogeneous beam equation. The consequences of a linearly varying pile distance ( ) were discussed in Section 2.1.4. By merely introducing ( ) rather than ( ) = , the analytical solutions became exceedingly complex. Therefore, we believe that even if a general solution was derived for non-uniform soil, it would most likely not be suited for a simplified, spreadsheet-friendly computational model. Even so, there are potential areas where the simplified models may be of value: 1. In the early stage of a design process or whenever approximate solution are sought. 2. Rough validation of more rigorous models. 3. In practical design situations, it is often interesting to observe trends rather than response values with respect to certain parameters. 4. Simplified methods are generally well-suited for educational purposes. Rigorous finite element methods are able to produce accurate results for a variety of structural and geotechincal problems. For better or worse, numerical models tend to directly embody and thus conceal some of the phenomena related to geodynamics and soil-structure-interaction, which are often important in the pursuit of mastering the above-mentioned fields of engineering. The derivation of the analytical model elucidates the majority of aspects related to structural dynamics, geodynamics, soil-pile interaction, computational mechanics and their interplay with each other.

Extension to the hybrid model
The hybrid model is constructed in Matlab [35] using beams, Winkler complex springs and especially formulated analogue interaction elements based on the attenuation functions given in Equations 15-18. The finite element stiffness and mass matrices are based on uniform, finite elements in accordance with the Euler-Bernoulli beam theory and the damping matrix is obtained as Rayleigh damping. [36] The Winkler springs and dashpots are modelled as SDOF-elements in accordance with Equations 1 and 3 and are directly incorporated as additional element matrices. Local element matrices are transformed to the global coordinate system trough a transformation matrix. In order to express the pile-soil-pile interaction, an analogue element is formulated by utilizing the plain strain attenuation function and the complex Winkler stiffness. The analogue element connects the horizontally aligned nodes in each beam, and the interaction is based on the plane-strain attenuation function and the geometrical considerations discussed in Sections 2.1.3 -2.1.7. See illustration in Figure 21. By considering transverse, rotational and axial degrees of freedom in each node, the complex stiffness matrix  . The analogue element stiffness matrix is complex since damping is directly incorporated into the matrix entries. The global stiffness matrix is thus also complex. Alternatively, the analogue element stiffness matrix may be split into separate stiffness and damping matrices, both of which are real. for example, ,1 may be split into stiffness and damping contributions, that is, The analogue element matrices are then treated in the same manner as the local finite element matrices and assembled into their corresponding global matrices, where both and only contain real entries. Figure 22 presents the horizontal, dynamic group stiffness for a two-by-one pile group with batter angles of 10 and 15 degrees and ∕ -ratio equal to 200. The hybrid and analytical model match well. Some discrepancy is observed for close pile spacing. The results indicate that the simplifications regrading (1) the neglection of receiver pile damping when deriving the interaction factor, (2) elimination of negligible parts of the non-homogeneous differential equation solution and (3) the assumption of a constant attenuation function, are in fact reasonable. Nonetheless, there are two obvious advantages to the hybrid model compared to the analytical model. Firstly, the hybrid model allows for arbitrary pile length and boundary conditions (fixed-head, hinged-head, frictional and end-bearing) whereas the analytical model demands an unique solution for every different combination of boundary condition. Secondly, the hybrid model may easily be extended to include any of type superstructure, and thus allows for direct solutions of earthquake loading problems. It should be emphasized that the hybrid model utilizes frequency dependent expressions, such that any arbitrary loading must be divided into harmonic components. The corresponding solutions must then be superimposed.

CONCLUSION
2. The analytical model is able to represent the trends associated with pile distance and frequency rather well. 3. The presented hybrid model combines the adaptability of a finite element model with the simplicity of an analytical model by replacing the soil elements with (1) Winkler elements and (2) especially developed analogue interaction elements to account for pile-soil-pile interaction based on simple attenuation functions. For the examined boundary conditions and batter angles, the analytical model matches the hybrid model fairly well.

A C K N O W L E D G M E N T S
The first author would like to express gratitude to Sweco Norway and The Research Council of Norway for financial support.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

C O N F L I C T O F I N T E R E S T
The authors declare no potential conflict of interests.