Flexible Pile Group Interaction Factors under Arbitrary Lateral Loading in Sand

Marine and harbor structures, wind turbines, bridges, offshore platforms, industrial chimneys, retaining structures etc. can be subjected to significant lateral loads from various sources. Appropriate assessment of the foundations capacity of these structures is thus necessary, especially when these structures are supported by pile groups. The pile group interaction effects under lateral loading have been investigated intensively in past decades, and the most of the conducted studies have considered lateral loading that acts along one of the two orthogonal directions, parallel to the edge of pile group. However, because of the stochastic nature of its source, the horizontal loading on the pile group may have arbitrary direction. The number of studies dealing with the pile groups under arbitrary loading is very limited. The aim of this paper is to investigate the influence of the arbitrary lateral loading on the pile group response, in order to improve (extend) the current design approach for laterally loaded pile groups. Free head, flexible bored piles in sand were analyzed through the extensive numerical study. The main hypothesis of the research is that some critical pile group configurations, loading directions, and soil conditions exist, which can lead to the unsafe structural design. Critical pile positions inside the commonly used pile group configurations are identified with respect to loading directions. The influence of different soil conditions was discussed.


Introduction
In situations where surface soil layers have low bearing capacity, pile foundations are the common foundation solution. Pile foundations are usually designed as closely spaced, square, or rectangular pile groups. Beside the primary function to transfer the vertical loads to the stiff soil layers, pile foundations can be significantly loaded with horizontal (lateral) loads. These loads may originate from different sources, such as: wave, current and ice action, ship impacts, wind pressure, earthquake, earth pressures, traffic acceleration, braking forces, soil displacements etc. The examples of laterally loaded engineering structures are: marine and harbor structures, wind turbines, bridges, offshore platforms, industrial chimneys, retaining structures, etc. The magnitude of lateral load is usually 10-15% of vertical load (and up to 30% in offshore structures) [1]. This fact makes the problem of the laterally loaded pile groups very challenging in marine engineering applications.
When the pile group is laterally loaded, stress-strain fields of adjacent piles overlap, followed with the separation ("gapping") between the piles and the soil at the back of the piles. The influence of leading (front) pile row on lateral response of trailing (rear) pile rows is called "shadowing" because the rear rows are in the shadow of the front pile row (Figure 1a). The problem of interaction between the piles inside the closely spaced laterally loaded pile group has been recognized by the between the piles inside the closely spaced laterally loaded pile group has been recognized by the scientific community. Ideally, the bearing capacity of the pile group shall be equal to the sum of the bearing capacities of the individual piles. However, due to the soil-structure interaction effects, loaddisplacement behavior of a pile inside the group is different from the behavior of the equivalent single pile (Figure 1b). In addition, the maximum pile bending moment in the pile group will be larger than for an equivalent single pile. When the pile spacing is increased, the effects of overlapping between the resisting zones become less significant. Engineering problems that can arise due to the inappropriate assessment of the pile lateral load can be very serious. Failure in such cases can lead to significant damage or even collapse of the entire superstructure. In the cases of bridges or other structures supported by pile foundations, only a few centimeters of lateral displacements can cause the significant stress development [2]. According to Eurocode 7 provisions [3], lateral pile displacements are usually limited to 3% of pile diameter, or to max. 2 cm. Despite the fact that the current design practice for pile foundations is based on the capacity-based design, many authors [2,[4][5][6][7][8][9] pointed out that the opposite, displacement-based design is more appropriate.
The pile group interaction effects under lateral loading have been investigated intensively in past decades, both numerically and experimentally. Many numerical solutions with different level of complexity have been proposed: closed form and empirical solutions [10][11][12], limit equilibrium methods [13,14], strain wedge method [15][16][17][18][19], and p-y curve method [20][21][22][23][24][25][26]. The development of fast computers led to rapid development of continuum-based methods, which allow for the full discretization of both soil and the structure, combined with the use of advanced constitutive models. Continuum-based representation of both the soil and the piles can provide a more realistic analysis, with the use of material properties with a clear physical meaning. On the contrary, continuum-based methods require more engineering effort for the numerical model preparation, as well as the higher computational hardware requirements. Continuum-based studies on the laterally loaded pile groups are mainly based on the following numerical methods: finite element method (FEM) [27][28][29][30][31], finite difference method (FDM) [2,8,[32][33][34] and boundary element method (BEM) [35][36][37].
Beside numerical studies, the pile group behavior under lateral loading has been investigated using different experimental methods. Full scale tests can eliminate most of the uncertainties of the considered problem. However, due to the high costs, only a limited number of full scale pile load test results are available in the literature. As an alternative, experiments with small scale models are also reported. The conducted experimental studies of the laterally loaded pile groups are summarized in several studies (e.g., Fayyazi et al. [38], Lemnitzer et al. [39], Ashour and Ardalan [40]). The results are mainly presented as the p-multipliers for various pile groups. Recently, improvement of CT (computer tomography) and PIV (particle image velocimetry) technologies led to the experimental studies [41][42][43][44][45] of soil deformation patterns under lateral loading, with the main conclusions Engineering problems that can arise due to the inappropriate assessment of the pile lateral load can be very serious. Failure in such cases can lead to significant damage or even collapse of the entire superstructure. In the cases of bridges or other structures supported by pile foundations, only a few centimeters of lateral displacements can cause the significant stress development [2]. According to Eurocode 7 provisions [3], lateral pile displacements are usually limited to 3% of pile diameter, or to max. 2 cm. Despite the fact that the current design practice for pile foundations is based on the capacity-based design, many authors [2,[4][5][6][7][8][9] pointed out that the opposite, displacement-based design is more appropriate.
The pile group interaction effects under lateral loading have been investigated intensively in past decades, both numerically and experimentally. Many numerical solutions with different level of complexity have been proposed: closed form and empirical solutions [10][11][12], limit equilibrium methods [13,14], strain wedge method [15][16][17][18][19], and p-y curve method [20][21][22][23][24][25][26]. The development of fast computers led to rapid development of continuum-based methods, which allow for the full discretization of both soil and the structure, combined with the use of advanced constitutive models. Continuum-based representation of both the soil and the piles can provide a more realistic analysis, with the use of material properties with a clear physical meaning. On the contrary, continuum-based methods require more engineering effort for the numerical model preparation, as well as the higher computational hardware requirements. Continuum-based studies on the laterally loaded pile groups are mainly based on the following numerical methods: finite element method (FEM) [27][28][29][30][31], finite difference method (FDM) [2,8,[32][33][34] and boundary element method (BEM) [35][36][37].
Beside numerical studies, the pile group behavior under lateral loading has been investigated using different experimental methods. Full scale tests can eliminate most of the uncertainties of the considered problem. However, due to the high costs, only a limited number of full scale pile load test results are available in the literature. As an alternative, experiments with small scale models are also reported. The conducted experimental studies of the laterally loaded pile groups are summarized in several studies (e.g., Fayyazi et al. [38], Lemnitzer et al. [39], Ashour and Ardalan [40]). The results are mainly presented as the p-multipliers for various pile groups. Recently, improvement of CT (computer tomography) and PIV (particle image velocimetry) technologies led to the experimental studies [41][42][43][44][45] of soil deformation patterns under lateral loading, with the main conclusions regarding the three-dimensional (almost conical) shape of the soil deformation zone, as well as the progressive soil failure.
The magnitude of the pile group interaction can be described by comparing the pile group response with the response of equivalent single pile for the same mean loading. One of the possible representations is the concept of pile interaction factors, which is still in use in everyday engineering practice in Germany, through the recommendations of DIN/EN codes [46,47] and working group "piles" (EAP) [48]. The pile interaction factor α i for the pile i inside the pile group is defined using the following expression: where H i is the lateral force on pile i and H SP is the lateral force for on equivalent single pile, for the same displacement level. According to the studies by Klüber [49] and Kotthaus [50], pile interaction factors are mainly dependent on both the pile spacing and pile position inside the pile group.
Interaction factors α i are calculated using the partial interaction factors α L (in the loading direction), α QA (perpendicular to loading direction, outer piles) and α QZ (perpendicular to loading direction, inner piles), using following equations: The determination of pile interaction factors is illustrated in Figure 2. regarding the three-dimensional (almost conical) shape of the soil deformation zone, as well as the progressive soil failure. The magnitude of the pile group interaction can be described by comparing the pile group response with the response of equivalent single pile for the same mean loading. One of the possible representations is the concept of pile interaction factors, which is still in use in everyday engineering practice in Germany, through the recommendations of DIN/EN codes [46,47] and working group "piles" (EAP) [48]. The pile interaction factor αi for the pile i inside the pile group is defined using the following expression: where Hi is the lateral force on pile i and HSP is the lateral force for on equivalent single pile, for the same displacement level. According to the studies by Klüber [49] and Kotthaus [50], pile interaction factors are mainly dependent on both the pile spacing and pile position inside the pile group. Interaction factors αi are calculated using the partial interaction factors αL (in the loading direction), αQA (perpendicular to loading direction, outer piles) and αQZ (perpendicular to loading direction, inner piles), using following equations: The determination of pile interaction factors is illustrated in Figure 2.  [50]). sL and sQ denote the pile center-to-center spacings in the direction of the lateral loading HG, and perpendicular to the direction of the lateral loading, respectively. Four different pile types, based on the pile position inside the pile group, are distinguished and associated with four different colors.
Proposed values for the partial interaction factors αL, αQA, and αQZ are given in Figure 3. Based on the recommendation of EAP [48], the pile interaction factors αi are used for the assessment of the modulus of subgrade reaction and pile bending moments for each pile inside the pile group.
(a) pile spacings (b) pile spacing perpendicular  [50]). s L and s Q denote the pile center-to-center spacings in the direction of the lateral loading H G , and perpendicular to the direction of the lateral loading, respectively. Four different pile types, based on the pile position inside the pile group, are distinguished and associated with four different colors.
Proposed values for the partial interaction factors α L , α QA , and α QZ are given in Figure 3. Based on the recommendation of EAP [48], the pile interaction factors α i are used for the assessment of the modulus of subgrade reaction and pile bending moments for each pile inside the pile group.
regarding the three-dimensional (almost conical) shape of the soil deformation zone, as well as the progressive soil failure.
The magnitude of the pile group interaction can be described by comparing the pile group response with the response of equivalent single pile for the same mean loading. One of the possible representations is the concept of pile interaction factors, which is still in use in everyday engineering practice in Germany, through the recommendations of DIN/EN codes [46,47] and working group "piles" (EAP) [48]. The pile interaction factor αi for the pile i inside the pile group is defined using the following expression: where Hi is the lateral force on pile i and HSP is the lateral force for on equivalent single pile, for the same displacement level. According to the studies by Klüber [49] and Kotthaus [50], pile interaction factors are mainly dependent on both the pile spacing and pile position inside the pile group. Interaction factors αi are calculated using the partial interaction factors αL (in the loading direction), αQA (perpendicular to loading direction, outer piles) and αQZ (perpendicular to loading direction, inner piles), using following equations: The determination of pile interaction factors is illustrated in Figure 2. Proposed values for the partial interaction factors αL, αQA, and αQZ are given in Figure 3. Based on the recommendation of EAP [48], the pile interaction factors αi are used for the assessment of the modulus of subgrade reaction and pile bending moments for each pile inside the pile group.
(a) pile spacings (b) pile spacing perpendicular Total pile group efficiency under lateral loading can be described using the following non-dimensional parameter (n-total number of piles) GW: The pile group efficiency GW is equal to the mean value of all interaction factors α i . Lateral force distribution between the piles can also be described using the load proportions LP i : So far, most of the studies on pile group interaction effects have considered lateral loading acting along one of the two orthogonal directions, parallel to the edge of pile group (namely X and Y directions, Figure 1a). However, because of the stochastic nature of its source, the horizontal loading on the pile group may have arbitrary direction. The number of studies dealing with the pile groups under arbitrary loading is very limited. Randolph [51] defined the influence of the loading direction via simple algebraic expressions for two laterally loaded flexible piles, both for elastic homogeneous and Gibson soil. Ochoa and O'Neill [52], based on the experimental results in submerged sand, have developed the interaction factors dependent on the angle between the load vector and the line that connects the pile heads. Fan and Long [53] analyzed the interaction between the piles under various loading directions and derived the modulus reduction factors to account for the interaction effects at the ultimate limit state. Su and Yan [54] formulated the multidirectional p-y model for sands that was incorporated into FEM and validated through the simulations of piles under unidirectional and multidirectional lateral loading. Mayoral et al. [55] formulated the p-y curves for piles under multidirectional loading in soft clays. The independent application of the common p-y curves in two orthogonal directions was found to be impossible. Georgiadis et al. [56] analyzed the influence of the arbitrary lateral loading on the ultimate resistance of the two rigid piles using FEM.
In their recent paper, Su and Zhou [57] presented the experimental study of 2 × 2 square and rectangular pile groups under static lateral loading in different directions. The results of this study have shown that the loading direction has great influence on the response of the pile group, mainly on the redistribution of the lateral load between the piles, and on the total group bearing capacity. The pile groups at medium pile spacing were found to be more sensitive to the variation in loading direction. In addition, some pile load distributions have increased when the loading direction was different than the usual 0/90 • case, which leads to the hypothesis that some critical (worst-case) loading cases exist, which can lead to the minimum foundation capacity. Such findings of this study, as well as the fact that, up to the authors' knowledge, arbitrary loading case was studied in a very limited number of research papers, mainly motivated further research presented in this paper.
Therefore, the main aim of this paper is to investigate the influence of arbitrary lateral loading on the pile group response, in order to improve (extend) the current design approach for laterally loaded pile groups. The main hypothesis of the research is that some critical pile group configurations, loading directions, and soil conditions exist, which can lead to unsafe design. This hypothesis is analyzed through the broad numerical study, and critical pile positions inside the commonly used pile group configurations are identified with respect to loading directions. The influence of different soil conditions, mainly sandy soil, is discussed.

Research Methodology and Scope
The methodology used in this study follows the concept of "numerical experiment", where a numerical simulation is used to mimic the real experiment. Compared to the real experiments, numerical simulations allow for the large number of analyses to be done faster and cheaper. Full scale experiments are expensive and therefore not feasible solution for the larger studies. First, the numerical model of laterally loaded pile group was generated and validated, and then the broad parametric study was conducted in order to assess the influence of various parameters on the pile group interaction.
Interaction inside the pile group is influenced by various factors, such as: pile spacing and arrangement inside the pile group, total number of piles, pile stiffness, length and diameter, pile installation method, soil conditions, loading level, pile head conditions, pile-soil slippage and separation, and presence of axial loading. A significant number of research papers [58][59][60][61][62][63][64][65] points out and ranks the most important factors that influence this problem. "Apriori" sensitivity analysis of different problem parameters, through the broad literature survey, was done in order to select the most important parameters for further analysis.
Regarding the scope of the present study, bored piles are analyzed, as more expected in common engineering practice. These piles can also sustain larger vertical loading, so the horizontal loading is expected to be high as well. According to Van Impe [66], "bored and CFA piles account for 50% of the world pile market, while the remaining is mainly covered by driven (42%) and screw piles". Regarding the pile length, long, flexible piles are considered, as the more common case [51]. Pile spacing and configuration are identified as the governing factors for the pile group interaction by different authors [2,67], while pile length and diameter are identified as the less important [32,51].
In this study, only free head piles were considered. According to Mokwa and Duncan [68], boundary conditions at the pile head are somewhere in between the fixed and free head pile. The same authors pointed out that the pure fixed head conditions are hard to achieve in reality, even when the pile cap is very stiff. Despite the different deflection response of fixed and free head piles, pile load redistribution can be considered independent on pile head conditions, especially within the working load conditions. This statement is supported by various authors [21,32,69].

Numerical Model
The FEM is today considered as the most reliable and widely used numerical method for engineering analysis of complex foundation systems. Regarding the laterally loaded pile group, FEM allows for detailed modeling of all important model components: pile geometry, soil continuity and nonlinear behavior, and especially the pile-soil interaction through the slippage and gapping. According to Dodds and Martin [23], "modeling lateral behavior in any way other than with three-dimensional models using nonlinear soil models and interface elements must constitute a compromise." Eurocode 7 [3] strongly suggests the use of soil-structure interaction in the numerical analysis, especially in the case of complex foundations, such as laterally loaded pile foundations or piled-raft foundations.
Within this study, the response of the pile groups under arbitrary static lateral loading was computed using FEM code PLAXIS 3D (Anniversary Edition) [70]. PLAXIS 3D is a world-wide used code for the stress-strain, stability, and groundwater flow analysis in geotechnical engineering, that supports an easy graphical input of the models with complex geometry, as well as the illustrative presentation of the results. It also features various constitutive models for soil and rock, with different levels of complexity.

Model Components
Numerical prediction of the laterally loaded pile group behavior is a challenging task that requires realistic modeling of the deformation behavior of the soil, pile, and the contact between the soil and the pile. As follows, model details for different model components are briefly presented.

Pile
Two pile modeling techniques are mainly used in PLAXIS 3D: full 3D (solid) pile model, and recently implemented embedded beam model. The comparison of two modeling techniques was studied by Marjanović et al. [71] on an idealized example of 2 × 2 pile group in loose and dense sand, with conclusions that, in the case of the weaker pile-soil contact, pile group behavior cannot be resembled using the embedded beam model in PLAXIS 3D. Therefore, the full 3D pile model was adopted for this study. Pile volumes are modeled using 10-node tetrahedral elements (nodes in the corners and the middles of the edges), with three translational degrees of freedom (u x , u y , u z ) per node. This type of element provides a second-order (quadratic) interpolation of nodal displacement components.
Pile behavior is assumed as linear elastic. Bowles [72] pointed out the possibility of concrete cracking as the problem in realistic modeling of the laterally loaded pile. Comodromos et al. [73] analyzed the influence of concrete cracking on the response of laterally loaded piles, and concluded that, at the working load levels, fully elastic assumption of the pile behavior will not lead to the significant change in the interaction effects and stress redistribution. Therefore, the pile constitutive behavior is defined using the linear elastic model, based on the generalized Hooke's law.
The extraction of the pile group results from the 3D FEM simulation is not a simple task. This is mainly related to the determination of pile section forces (shear forces and bending moments). Theoretically, this task is done by integration of component stresses over the pile cross section. Despite the fact that the calculation of the section forces from volume elements is implemented in PLAXIS 3D, this step cannot be done automatically, but requires that the user performs the whole procedure manually for each pile volume. In order to simplify the calculation of pile section forces and to automatize the post-processing as well, the "dummy" beams are added into the numerical model. Such a modeling concept is common when the laterally loaded piles are modeled using the 3D elements [74,75]. The elastic beam finite elements with very small bending stiffness (10 6 times smaller than the pile bending stiffness EI) are inserted along the pile axes. Because the "dummy" beam stiffness is very small, system stiffness matrix remains almost unchanged. "Dummy" beams are enforced to deform together with the pile axis, and the pile section forces are then easily computed by simple multiplication of "dummy" beam section forces with 10 6 .

Pile-Soil Interface
In reality, the contact between the piles and the surrounding soil is not fully rigid, but weaker, so there are small relative displacements along the pile-soil interface. In general, the pile-soil interface strength depends on the type of pile material (wood, steel, concrete, etc.) and pile installation effects (bored, driven, jacked-in). Because this study considers only the bored piles, installation effects and loading history were neglected.
The pile-soil contact was modeled using thin 2D interface elements. These elements are different from the regular finite elements: they have pairs of nodes instead of single nodes, with the distance between the two nodes of a node pair equal to zero. Each node has three translational degrees of freedom (u x , u y , u z ). As the result, these elements allow for differential displacements between the node pairs to simulate both slipping and gapping on the pile-soil contact [70].
The constitutive behavior of the pile-soil interface is defined by the elastic-perfectly plastic Mohr-Coulomb (MC) model with a non-associated flow rule and zero tension cut-off criterion (when tension develops, a gap between the pile and the surrounding soil is generated). The shear strength parameters of the pile-soil interface (c inter , ϕ inter ) are related to the strength reduction factor R inter , which reduces the shear strength parameters (c and ϕ ) of the MC model, according to the following expressions: tan ϕ inter = R inter tan ϕ .
In general, there is lack of experimental data for real pile-soil interface parameters. Suitable values of R inter are recommended in the literature [70] for different soil types, and the usual value is around 0.5.

Soil
Soil constitutive behavior is modeled using the Hardening Soil (HS) model (Schanz and Vermeer [76]). As in the case of the MC model, the limit stress states of the HS model are defined using the parameters for the MC failure criterion (cohesion c , angle of internal friction ϕ , and dilation angle ψ). However, the soil stiffness is described more accurately, using advanced stiffness parameters: the triaxial loading stiffness E 50 , the triaxial unloading stiffness E ur , and the oedometer loading stiffness E oed . These stiffnesses are not constants, but are dependent on the loading level (principal stress state). The magnitude of the stress level dependency is governed using the parameter m. The values of these three stiffness parameters are not mutually independent-usually the following relationships are used: Opposite to the MC model, the yield surface of the HS model is not fixed in the principal stress space, but it can expand due to plastic straining. Two types of hardening are included in the model: shear hardening (due to primary deviatoric loading) and compression hardening (due to oedometer and isotropic loading).
Due to the fact that the pile group response under lateral loading is mainly governed by the properties of the top soil layer (near-surface soil conditions) [2], the soil was modeled as homogeneous and isotropic. Because the static loading is assumed, drained soil conditions are assumed in the model. The soil volumes are also modeled using 10-node tetrahedral elements.

Model Geometry, Discretization, and Boundary Conditions
Model side domain limits were chosen far enough from the center of the pile group, in order to avoid the boundary effects on the model response. Sensitivity analyses were done to determine the optimum size of the model. Because the model geometry will change for different pile group configurations, it was adopted that the distance between the model domain sides and the edges of the pile group remain the same. The spacing between the pile bottom (tip) and the bottom model boundary was chosen to be relatively small, because the lateral pile group response is mainly governed by the top soil layer in the case of flexible piles, as well as with the active pile length.
PLAXIS 3D does not allow for very precise control of the model discretization. However, it allows the user to define "refinement zones" inside the model volume, where the mesh discretization can be finer. Around the pile group, the rectangular prismatic refinement zone was adopted. Trial analyses were executed on models with different discretization setups, in order to optimize the accuracy and computational cost. The most of the model geometry parameters have been designed as the multiplicators of pile diameter D, so the whole geometry can be parameterized as the function of pile diameter D.
The displacements on the model boundaries were fully fixed in all directions at the bottom of the model, while the displacements on the side planes were limited to vertical direction. The final model geometry and the mesh refinement zone are presented in Figures 4 and 5.

Calculation Stages
Numerical simulations within the study were done as the staged (phase) analysis, and consisted of the following stages:

1.
Initial (K 0 ) stage-initial stress field in the soil is established. Horizontal stress state is calculated using the Jaky's formula [77]: Construction stage-soil volume is replaced with 3D finite elements that represent the piles (wished-in-place concept). Interface elements are also activated in this stage 3.
Prescribed displacements (4 increments)-the prescribed displacements up to 0.1D are applied at the pile top, in desired loading direction. The prescribed displacements simulate the displacement control test under static loading conditions.
2. Construction stage-soil volume is replaced with 3D finite elements that represent the piles (wished-in-place concept). Interface elements are also activated in this stage 3. Prescribed displacements (4 increments)-the prescribed displacements up to 0.1D are applied at the pile top, in desired loading direction. The prescribed displacements simulate the displacement control test under static loading conditions.   2. Construction stage-soil volume is replaced with 3D finite elements that represent the piles (wished-in-place concept). Interface elements are also activated in this stage 3. Prescribed displacements (4 increments)-the prescribed displacements up to 0.1D are applied at the pile top, in desired loading direction. The prescribed displacements simulate the displacement control test under static loading conditions.

Post-Processing of the Results
Beside the complexity of described numerical model and numerical simulation itself, data management after calculations was an equally challenging task, because all raw data had to be reordered to interpret the pile group results. PLAXIS 3D allows the extraction of all results in Euclidean XY space. However, for the arbitrary loading case, the resultant maximum displacement, shear forces, and bending moments must be recalculated, using simple vector algebra ( Figure 6):

Post-Processing of the Results
Beside the complexity of described numerical model and numerical simulation itself, data management after calculations was an equally challenging task, because all raw data had to be reordered to interpret the pile group results. PLAXIS 3D allows the extraction of all results in Euclidean XY space. However, for the arbitrary loading case, the resultant maximum displacement, shear forces, and bending moments must be recalculated, using simple vector algebra ( Figure 6): The most important step in the post-processing of the pile group results is the precise evaluation of the pile shear forces. However, the shear forces extracted directly from the "dummy" beams in PLAXIS 3D are slightly unrealistic, which is the issue recognized and analyzed by Tedesco [74], who concluded that such behavior could be associated with the PLAXIS 3D beam elements, that compute the shear forces using the bending moment derivative along the pile length. In order to properly evaluate the pile shear forces, first, the pile bending moments were approximated using the B-spline approximations, with 10 interior knots and 5th order spline interpolation. Other fitting techniques, such as polynomials and cubic splines, were also considered. Then, the shear forces profile was computed by differentiation of the fitted pile bending moments line along the pile, and then the pile shear force at the ground level was used for the calculation of the interaction factors. The shape of the characteristic diagrams obtained in the numerical analysis is given in Figure 7.

Automatization of the Calculation Procedures
Parametric study done within this research includes multiple FEM simulations, with different combinations of input data. In such case, one of the critical "bottlenecks" are the repetitive pre-/postprocessing activities, and this part of the analysis is usually prone to errors. Despite the fact that the PLAXIS 3D is very user friendly, the preparation and (especially) the post-processing of multiple numerical models of the pile group is time-consuming. Scripting is a popular approach for the automatization of numerical analysis, and such features are nowadays implemented in PLAXIS 3D and many other FEM packages, mainly using the Python programming language. PLAXIS 3D input The most important step in the post-processing of the pile group results is the precise evaluation of the pile shear forces. However, the shear forces extracted directly from the "dummy" beams in PLAXIS 3D are slightly unrealistic, which is the issue recognized and analyzed by Tedesco [74], who concluded that such behavior could be associated with the PLAXIS 3D beam elements, that compute the shear forces using the bending moment derivative along the pile length. In order to properly evaluate the pile shear forces, first, the pile bending moments were approximated using the B-spline approximations, with 10 interior knots and 5th order spline interpolation. Other fitting techniques, such as polynomials and cubic splines, were also considered. Then, the shear forces profile was computed by differentiation of the fitted pile bending moments line along the pile, and then the pile shear force at the ground level was used for the calculation of the interaction factors. The shape of the characteristic diagrams obtained in the numerical analysis is given in Figure 7.

Post-Processing of the Results
Beside the complexity of described numerical model and numerical simulation itself, data management after calculations was an equally challenging task, because all raw data had to be reordered to interpret the pile group results. PLAXIS 3D allows the extraction of all results in Euclidean XY space. However, for the arbitrary loading case, the resultant maximum displacement, shear forces, and bending moments must be recalculated, using simple vector algebra ( Figure 6): The most important step in the post-processing of the pile group results is the precise evaluation of the pile shear forces. However, the shear forces extracted directly from the "dummy" beams in PLAXIS 3D are slightly unrealistic, which is the issue recognized and analyzed by Tedesco [74], who concluded that such behavior could be associated with the PLAXIS 3D beam elements, that compute the shear forces using the bending moment derivative along the pile length. In order to properly evaluate the pile shear forces, first, the pile bending moments were approximated using the B-spline approximations, with 10 interior knots and 5th order spline interpolation. Other fitting techniques, such as polynomials and cubic splines, were also considered. Then, the shear forces profile was computed by differentiation of the fitted pile bending moments line along the pile, and then the pile shear force at the ground level was used for the calculation of the interaction factors. The shape of the characteristic diagrams obtained in the numerical analysis is given in Figure 7.

Automatization of the Calculation Procedures
Parametric study done within this research includes multiple FEM simulations, with different combinations of input data. In such case, one of the critical "bottlenecks" are the repetitive pre-/postprocessing activities, and this part of the analysis is usually prone to errors. Despite the fact that the PLAXIS 3D is very user friendly, the preparation and (especially) the post-processing of multiple numerical models of the pile group is time-consuming. Scripting is a popular approach for the automatization of numerical analysis, and such features are nowadays implemented in PLAXIS 3D and many other FEM packages, mainly using the Python programming language. PLAXIS 3D input

Automatization of the Calculation Procedures
Parametric study done within this research includes multiple FEM simulations, with different combinations of input data. In such case, one of the critical "bottlenecks" are the repetitive pre-/post-processing activities, and this part of the analysis is usually prone to errors. Despite the fact that the PLAXIS 3D is very user friendly, the preparation and (especially) the post-processing of multiple numerical models of the pile group is time-consuming. Scripting is a popular approach for the automatization of numerical analysis, and such features are nowadays implemented in PLAXIS 3D and many other FEM packages, mainly using the Python programming language. PLAXIS 3D input and output modules can be fully controlled through the user defined Python scripts, instead the GUI (graphical user interface).
In order to perform the most of the calculation procedure in this study automatically, the broad development and implementation of original Python scripts was done throughout every numerical analysis step. In this manner, the PLAXIS 3D becomes a "problem calculator", while the simulation parameters and data analysis are handled by the separate computer programs. Such approach cannot increase the single simulation calculation time, neither prevent the conceptual errors in the numerical model, but data handling in all simulation stages becomes substantially faster.
Beside scripting, PLAXIS 3D allows for the use of multiple computers through the local area network. Due to availability of multiple licenses of PLAXIS 3D, the set of computer programs for multiple simulation running on different computers was developed. The parametric study was then done on (in average) 4-7 computers at the same time, depending on the current availability. The program was tested on up to 30 computers, and the only limit is the number of available PLAXIS 3D licenses. It is important to note that single computer calculates only one simulation at a time, but as soon the simulation is finished, the next simulation is started automatically. This process lasts as long as the simulation queue contains at least one simulation, and the total computation time is decreased with the number of used computers, assuming the same simulation complexity and the same computer hardware. This original solution has provided the optimum use of relatively expensive resources (code licenses and computer power).

Model Validation
Numerical model was validated based on the back-calculation of the small-scale pile group test from the literature. Trial and error analysis was used to match the experimental results with the model response.

Experimental Results
A well-documented centrifuge test at 50 g, done by Kotthaus [50], was used for model validation. Both single pile and the three pile rows at 3D spacing were back-calculated. The experimental setup was prepared by the sand rainfall method, to achieve the desired relative sand density. Dry sand was used for the centrifuge test. Experimental piles were fabricated as hollow aluminum tubes. The outside pile diameter was 30 mm, and the pipe wall thickness was 2 mm. Model dimensions were chosen to mimic the real structure, made of solid concrete piles. The model and prototype dimensions are given in Table 1. Prototype model dimensions are determined according to scaling laws [78]. Free head pile top conditions were provided using the screws that allowed for the rotation of the pile head, but enforced the same pile top displacements (displacements control test). In order to increase the pile soil friction inside the dense sand, soil particles were glued onto the pile shaft. The point of load application was above the ground level, because of the experimental setup.

Validation Results
For the accurate numerical simulation of geotechnical problem, the first step is selection and proper calibration of the soil constitutive model. The available experimental soil testing results for the dense sand (C U = 2.08, D 10 = 0.12 mm, D 50 = 0.23 mm, e max = 0.914, e min = 0.583) used in the centrifuge experiment by Kotthaus [50] were used for the calibration of the hardening soil model in PLAXIS 3D. The effective angle of internal friction for the dense sand was initially assumed as φ' = 41 • , and oedometer modulus as E oed = 28 MPa, based on the results of the sand laboratory tests given in [50,79]. However, slight trial and error adjustments of initially assumed values had to be done in order to better match the measured pile group response. This alteration is justified by the fact that the stress states around the laterally loaded are not fully matched by conventional laboratory tests. Finally adopted constitutive model parameters are given in Table 2. The measured vs. computed load-displacement response is given in Figure 8. Despite slight discrepancies, presented results of the model validation show a satisfactory match with the experimental results and the overall performance of the numerical model is considered to be acceptably accurate for the parametric study of the considered problem.

Validation Results
For the accurate numerical simulation of geotechnical problem, the first step is selection and proper calibration of the soil constitutive model. The available experimental soil testing results for the dense sand (CU = 2.08, D10 = 0.12 mm, D50 = 0.23 mm, emax = 0.914, emin = 0.583) used in the centrifuge experiment by Kotthaus [50] were used for the calibration of the hardening soil model in PLAXIS 3D. The effective angle of internal friction for the dense sand was initially assumed as ϕ' = 41°, and oedometer modulus as Eoed = 28 MPa, based on the results of the sand laboratory tests given in [50,79]. However, slight trial and error adjustments of initially assumed values had to be done in order to better match the measured pile group response. This alteration is justified by the fact that the stress states around the laterally loaded are not fully matched by conventional laboratory tests. Finally adopted constitutive model parameters are given in Table 2. The measured vs. computed load-displacement response is given in Figure 8. Despite slight discrepancies, presented results of the model validation show a satisfactory match with the experimental results and the overall performance of the numerical model is considered to be acceptably accurate for the parametric study of the considered problem.

Parametric Study
The scope of this study was designed with the aim to assess the most anticipated practical situations with the reasonable complexity. The following major problem parameters, according to "apriori" sensitivity analysis, were selected for this study: soil type, pile group configuration, and pile spacing, as well as the loading direction, as the main parameter under investigation (Figure 9). The selection of study parameters is briefly elaborated in the following subsections.

Parametric Study
The scope of this study was designed with the aim to assess the most anticipated practical situations with the reasonable complexity. The following major problem parameters, according to "apriori" sensitivity analysis, were selected for this study: soil type, pile group configuration, and pile spacing, as well as the loading direction, as the main parameter under investigation (Figure 9). The selection of study parameters is briefly elaborated in the following subsections. Two "synthetic" soil types, loose and dense dry sand, were simulated. Brown et al. [26] noted that the shadowing effect was more pronounced in sand compared with stiff clay. Constitutive model parameters for dense sand was adopted to be close to the parameters of the validated numerical model, and loose sand parameters were chosen as significantly different, in order to compare the soil stiffening effects on the pile group response. The parameters of the investigated soils are summarized in Table 3.  Two "synthetic" soil types, loose and dense dry sand, were simulated. Brown et al. [26] noted that the shadowing effect was more pronounced in sand compared with stiff clay. Constitutive model parameters for dense sand was adopted to be close to the parameters of the validated numerical model, and loose sand parameters were chosen as significantly different, in order to compare the soil stiffening effects on the pile group response. The parameters of the investigated soils are summarized in Table 3. Pile groups with different pile arrangements were analyzed, as summarized in Table 4. The numerical model symmetry could not be applied. In order to effectively analyze and discuss the interaction effects, as well as to limit the calculation time, it was decided to limit the total number of piles inside the pile group up to 3 × 3 piles. From the practical standpoint, it is a more common case that the configuration of the pile group (NxM) varies (square, rectangular, circular, etc.), while the spacing between the piles is kept on the lowest possible level. If the pile spacing is higher than 5D, the costs of the pile cap will increase to a very high level [60]. In this study, pile spacing between 2D-5D were investigated.

Results and Discussion
The results of the parametric study are presented in the form of pile interaction factors and maximum bending moments for each pile inside the pile group, with respect to loading direction. The representative plots are given to illustrate the main results trends. All results are presented for the displacement levels 0.03D (y/D = 0.03) and 0.10D (y/D = 0.10).
All original computer codes, PLAXIS 3D AE numerical models, raw simulation results and results databases are available from the first author via e-mail upon request.

Pile Interaction Factors
Pile interaction factors for different pile group configurations are presented on Figures 10-13. Full lines denote the results for loose sand, while dashed lines denote the results for dense sand.

Pile Interaction Factors
Pile interaction factors for different pile group configurations are presented on Figures 10-13. Full lines denote the results for loose sand, while dashed lines denote the results for dense sand.            Pile interaction factors for the closely-spaced square 2 × 2 pile group (s x = 2D, s y = 2D) are presented in Figure 10. The significant influence of the displacement level is observed-the interaction factors are changing in a broader range at the higher displacement level. The influence of the soil type is also pronounced, but at the smaller level, compared to the pile position inside the pile group. Similar trends can be seen for the square 3 × 3 pile group interaction factors, given in Figure 13. The effects of symmetry are clearly observed: the black and green pile have the same interaction factors for the loading direction angle equal to zero. For the direction angle of 45 degrees, the green and blue pile have the same interaction factors, due to symmetry. In addition, the well-known behavior of front and rear piles can be observed: for the loading direction angle equal to zero, black and green (front) piles carry higher loading than red and blue (rear) piles. Pile interaction factors for the closely-spaced square 2 × 2 pile group (sx = 2D, sy = 2D) are presented in Figure 10. The significant influence of the displacement level is observed-the interaction factors are changing in a broader range at the higher displacement level. The influence of the soil type is also pronounced, but at the smaller level, compared to the pile position inside the pile group. Similar trends can be seen for the square 3 × 3 pile group interaction factors, given in Figure  13. The effects of symmetry are clearly observed: the black and green pile have the same interaction factors for the loading direction angle equal to zero. For the direction angle of 45 degrees, the green and blue pile have the same interaction factors, due to symmetry. In addition, the well-known behavior of front and rear piles can be observed: for the loading direction angle equal to zero, black and green (front) piles carry higher loading than red and blue (rear) piles.
Pile interaction factors for the closely-spaced rectangular 2 × 2 pile group (sx = 2D, sy = 3D) are presented in Figure 11. The similar shape of the plot is also observed for different pile spacing, with decreased interaction with increased pile spacing. As in the previous case (Figure 10), the significant influence of the displacement level and medium influence of soil type on pile group interaction is observed. The behavior of green and blue piles is significantly different, compared to the other two piles. The changing of interaction factors for these piles is more pronounced with the variation in loading direction. Therefore, these piles can be considered as sensitive to change in loading direction, which can lead to the unsafe design in the case of unpredictable change of loading direction.
Pile interaction factors for the closely-spaced rectangular 2 × 3 pile group (sx = 2D, sy = 2D) are presented in Figure 12. The interaction factors in this case follows the same trends as in the previous case. The green, purple, and black pile can be considered as more sensitive.
Finally, pile interaction factors for square 3 × 3 pile group (sx = 3D, sy = 3D) are presented in Figure 13. The trend lines, in general, follows the same pattern as in the case of the 2 × 2 square pile group.

Maximum Bending Moments
Maximum pile bending moments for different pile group configurations are presented in Figures 14-16. Full lines denote the results for loose sand, while dashed lines denote the results for dense sand. Pile interaction factors for the closely-spaced rectangular 2 × 2 pile group (s x = 2D, s y = 3D) are presented in Figure 11. The similar shape of the plot is also observed for different pile spacing, with decreased interaction with increased pile spacing. As in the previous case (Figure 10), the significant influence of the displacement level and medium influence of soil type on pile group interaction is observed. The behavior of green and blue piles is significantly different, compared to the other two piles. The changing of interaction factors for these piles is more pronounced with the variation in loading direction. Therefore, these piles can be considered as sensitive to change in loading direction, which can lead to the unsafe design in the case of unpredictable change of loading direction.
Pile interaction factors for the closely-spaced rectangular 2 × 3 pile group (s x = 2D, s y = 2D) are presented in Figure 12. The interaction factors in this case follows the same trends as in the previous case. The green, purple, and black pile can be considered as more sensitive.
Finally, pile interaction factors for square 3 × 3 pile group (s x = 3D, s y = 3D) are presented in Figure 13. The trend lines, in general, follows the same pattern as in the case of the 2 × 2 square pile group.

Maximum Bending Moments
Maximum pile bending moments for different pile group configurations are presented in Figures 14-16. Full lines denote the results for loose sand, while dashed lines denote the results for dense sand.
As shown in Figure 14, maximum pile bending moments are not significantly influenced with the loading direction. As expected, higher bending moments occur in dense sand (due to the displacement control test and larger soil resistance in dense sand). Beside the higher absolute bending moment, the displacement level does not affect the distribution of the bending moments between the piles. It can be concluded that the pile group configuration and pile position are the governing factors in this case.   As shown in Figure 14, maximum pile bending moments are not significantly influenced with the loading direction. As expected, higher bending moments occur in dense sand (due to the displacement control test and larger soil resistance in dense sand). Beside the higher absolute bending moment, the displacement level does not affect the distribution of the bending moments between the piles. It can be concluded that the pile group configuration and pile position are the governing factors in this case.   As shown in Figure 14, maximum pile bending moments are not significantly influenced with the loading direction. As expected, higher bending moments occur in dense sand (due to the displacement control test and larger soil resistance in dense sand). Beside the higher absolute bending moment, the displacement level does not affect the distribution of the bending moments between the piles. It can be concluded that the pile group configuration and pile position are the governing factors in this case.   As shown in Figure 14, maximum pile bending moments are not significantly influenced with the loading direction. As expected, higher bending moments occur in dense sand (due to the displacement control test and larger soil resistance in dense sand). Beside the higher absolute bending moment, the displacement level does not affect the distribution of the bending moments between the piles. It can be concluded that the pile group configuration and pile position are the governing factors in this case. In the case of the rectangular 2 × 3 pile group (s x = 2D, s y = 2D), maximum pile bending moments ( Figure 15) are distributed more unequally between the piles, compared with the square pile group at the same pile spacing. Again, higher bending moments occur in dense sand. In addition, the influence of the lateral loading direction on the maximum bending moments can be considered as less important, compared to the pile position inside the pile group. The same pile groups (2 × 3) with different spacing show similar trends, with less moment differences between the piles at larger pile spacing (as expected).
Maximum pile bending moments for the square 3 × 3 pile group (s x = 3D, s y = 3D) follows the similar trend as the smaller square pile group, but with significantly higher differences of the maximum moment between the piles. In general, the maximum bending moment between the piles shows larger discrepancies with the increased number of the piles in the pile group configuration.

Conclusions
A broad numerical experiment has been conducted to examine the influence of the arbitrary lateral loading direction on the interaction inside the pile group. Free head, flexible bored piles in dry sands were considered. The FEM 3D numerical model was designed and validated by back-calculation of the existing experimental results. The capabilities of PLAXIS 3D software to perform multiple numerical simulations have been tested in an automated scripting environment, using original computer programs. The common approach to the analysis of the laterally loaded pile group is expanded with the pile interaction factors for the arbitrary loading direction. The results of the study provide more clear insight into the pile group behavior under lateral loading. Based on the presented analyses, the following conclusions can be made:

1.
Modeling of the laterally loaded pile group using the full 3D finite element model provides reasonable results for the case of the bored piles.

2.
As expected, the level of interaction between the piles is higher at higher displacement levels. The influence of the displacement level on the pile group interaction effects should be investigated in more detail, especially at small displacement levels.

3.
The interaction level is decreased with increased pile spacing. By means of quantification, interaction factors are between 0.6-0.9 for working load levels, and 0.4-0.9 for high loading levels. 4.
The pile interaction factors are dependent on soil type, but soil conditions are less important factor, compared to pile group configuration, displacement level, and pile position inside the pile group.

5.
The influence of the loading direction on the maximum bending moments is relatively small, which is expected, given the central symmetry of the circular piles. The soil conditions significantly affect the bending moments. The displacement level does not affect the maximum bending moment distribution between the piles. Bending moment discrepancies are more pronounced with the increased number of piles. 6.
The force in some individual piles inside the pile group significantly changes with the change of loading direction. These piles can be considered as sensitive, because the change of the loading direction can be unpredictable. The critical positions of these piles inside the pile group have been identified for the considered pile group configurations. 7.
The concept of multiple numerical simulations, followed with original computer programs, could be easily extended to other numerical models and more used in everyday practice. Bearing in mind that the simulation sets can be run during the night, this approach leads to the optimum use of hardware and software resources.
The presented conclusions emphasize the importance of the considered problem and can lead to further research topics in various fields. Due to availability of modern experimental techniques, such as CT and PIV, the experimental analysis of soil deformation patterns under arbitrary loading would provide better insight in the problem and lead to more accurate numerical models. The identified concept of sensitive piles should be investigated in more detail, with the aim to eliminate the sensitive elements from the foundations design. The influence of the different water level in the soil could also be investigated in more detail, since this paper considered only dry sand. Finally, the speed of numerical simulations could be reduced with the alternative modeling approaches-through development of improved embedded beam model formulations (with advanced interface) in FEM, as well as the more advanced strain wedge model (with arbitrary wedge orientation).