Practical macro‐element for vertical and batter pile groups

This paper presents a novel macro‐element for vertical and batter pile groups. The numerical scheme is based on de‐coupled, single pile response. Each pile consists of two separate load‐displacement formulations (axial and transverse) that take a displacement increment as input and return a diagonal tangent stiffness value. The effect of rotation is implicitly incorporated in the transverse load‐displacement formulation. The presented macro‐element does not require pre‐defined failure surfaces or other parameters, and is therefore not restricted to a specific foundation configuration, soil profile or soil type. The macro‐element may be calibrated using any type of non‐linear pile‐soil model. The developed macro‐element is validated using rigorous numerical models for a variety of configurations.

cyclic loading, soil plasticity and uplift. The element was later applied on the seismic analysis of a bridge pier foundation. 23 Martin and Houlsby 24,25 presented a monotonic model for spudcan footings on clay based on experimental tests. A similar model was developed by Houlsby and Cassidy 26 for sand. Cassidy et al. 27 presented a monotonic macro-element for spudcan footings with six degrees-of-freedom. Houlsby et al. 28 derived a theoretical model based on Winkler springs for quasi-static loading. Einav and Cassidy 29 developed a model similar to Houlsby et al. 28 Salciarini and Tamagnini 30 formulated a model for shallow footings on sand which shared some similarities with the work of Nova and Montrasio. 17 Chatzigogos et al. 31 presented a bounding plasticity model considering soil inelasticity and non-linear uplift mechanism. Chatzigogos et al. 32 extended this model to include sliding. Figini et al. 33 expanded on the work of Chatzigogos et al. 31,32 and validated the results using cyclic and dynamic large-scale laboratory tests. Ibsen et al. 34 and Foglia et al. 35 investigated the response of bucket foundations on dry sand. Skau et al. 36 presented a macro-element for bucket foundations based on multi-surface plasticity. Tistel and Grimstad 37 presented a monotonic model for an anchor block foundation on sand. Millen et al. 38 presented a macro-element for shallow foundations based on the work of Chatzigogos et al. 31 and Figini et al. 33 .
In recent years, attempts have been made to develop macro-elements for deep foundations. Correia 39 and Correia and Pecker 40 presented a pile-head macro-element for monoshaft foundations. The formulation accounted for non-linear behaviour of pile, soil and separation effects. In addition, a rigorous calibration procedure was presented. Inspired by the work of Salciarini and Tamagnini, 30 Liu et al. 41,42 first developed a macro-element for single piles embedded in homogeneous sand, which they later extended to single batter piles. Page et al. 43 presented a macro-element model for mono-pile foundations based on multi-surface plasticity and verified it against large-scale model tests. 44 Pérez 45 presented a macro-element for vertical pile groups based on the work of Liu. 46 To the authors knowledge, there are no macro-elements developed for pile groups with vertical and batter piles that take into account the inelastic behaviour of both pile and soil. The objective of this paper is to formulate a practical macroelement with three degrees of freedom for vertical and batter pile groups. Since macro-elements are most suitable for practical engineering and design tools, the use of such elements should be more attractive than finite element modelling of the entire system. Consequently, certain criteria will be emphasized in the formulation. First, the macro-element does not have to capture all the features of a rigorous finite element model, but it needs to capture the trends intrinsic to pile group configuration and soil profile. Second, the formulation must be robust with respect to both pile group configuration and numerical implementation. Third, the calibration must be straight-forward and independent of pile and soil properties, soil profile and pile group configuration. To achieve robustness with respect to pile group configuration and straightforward calibration procedures, the numerical scheme presented in this study is based on de-coupled, single pile response. Each pile consists of two separate load-displacement formulations (axial and transverse) that take a displacement increment as input and return a diagonal tangent stiffness value. The effect of rotation is implicitly incorporated in the transverse load-displacement formulation. The global tangent stiffness matrix (which is passed to the global solution in a finite element code) is assembled on the basis of the single pile tangent stiffness values. The presented macro-element does not require pre-defined failure surfaces or other parameters, and is therefore not restricted to a specific foundation configuration, soil profile or soil type. It should be noted that even though the calibration in this paper is achieved using detailed finite element models, the macro-element may be calibrated using any type of non-linear pile-soil model. However, there are limitations that must be addressed. First, the macro-element is formulated for seismic design purposes where large displacements are expected. It is assumed that geometric damping (frequency-dependent) is negligible compared to material damping (frequency-independent). The constitutive model is therefore rate-independent. Second, pile-soil-pile interaction is neglected. This assumption admittedly decreases the accuracy of the proposed model when the piles are closely spaced. While it is possible to approximately include pile-soil-pile interaction, it should be noted that recent studies 47,48 have shown that pile-to-pile interaction is less significant when piles undergo large displacements in soft, inelastic soil. Third, the macro-element is formulated such that pile-cap rotation is resisted by vertical pile-head loads only. This assumption may also decrease the accuracy when the piles are closely spaced. Finally, axial and transverse response are de-coupled. It is assumed that the active length contributes to transverse resistance while the soil below the active length contributes to axial resistance. This assumptions strictly limits the formulation to long piles. The effect of axial load on bending capacity is not considered. This paper is divided in five sections. Section 2 briefly describes the finite element model used for validation. Section 3 presents the theoretical framework, numerical treatment and validation of the single pile macro-elements. Section 4 outlines the pile group tangent stiffness matrix assembly and demonstrates the pile group macro-element performance for numerous configurations. The conclusion is given in Section 5.
(A) Cross section of the concrete pile (B) Soil profile.

VALIDATION MODEL
The macro-element is validated by finite element models constructed in OpenSees MP 49 together with the pre-and postprocessing tool STKO. 50 The soil profile and cross section of the concrete pile are shown in Figure 1. The piles are modelled using displacement-based beam elements. The non-linear behaviour of reinforced concrete is represented using fibre sections. The confined and un-confined concrete is simulated using the uni-axial material models ConfinedConcrete01 51 and Concrete01, 52 respectively. The cylindrical strength is 45 and the initial elastic stiffness modulus is 36 . The reinforcement is simulated using the uni-axial material model Steel02. 53 The yield strength is 487 and the initial elastic stiffness modulus is 185 . The soil is divided in six layers along the pile length, where each layer has a height ℎ = 3 . The soil is modelled using eight-noded hexahedral elements with a single integration point to prevent locking behaviour. The adopted soil model Pressure Independent Multi-Yield is an elastic-plastic, soil model suited for clay. 54 The reader is referred to Cemalovic et al. 55 for a detailed description of the soil material model, soil profile and meshing strategies. The pile-soil interface is modelled using frictional contact elements based on the Mohr-Coulomb criterion, penalty constraints and an implicit-explicit solution scheme. 50,56 Multi-stage analyses are performed to capture the correct stress state prior to the quasi-static pile-head loading. The first stage is a gravity analysis of the soil domain. In the second stage, a new gravity analysis is performed where the soil corresponding to the pile geometry is removed. In the third stage, a new gravity analysis is performed which includes the piles and contact elements. In the fourth and final stage, the quasi-static pile-head load analysis is performed. Prior to each analysis, the displacement field is set equal to zero. The system is solved using the Krylov-Newton implicit scheme. 57

One-dimensional bounding plasticity
Several researchers have successfully implemented bounding surface plasticity models for shallow and deep foundations. [31][32][33]36,39,40,43 The macro-element presented in this paper is based on the bounding surface plasticity theory with radial mapping. 58,59 There are three major aspects that distinguish bounding surface plasticity from conventional rate-independent plasticity. First, the current load is limited by a bounding surface. Second, the current load surface is also the yield surface, which implies that a purely elastic domain does not exist (except for a infinitesimal range at initial loads or load reversals). Third, the evolution of internal variables is governed only by the distance from the current load to a image point on the bounding surface, which implies that there is no need for an explicit description of a hardening variable. Since the macro-element in this study considers un-coupled single pile response, we will henceforth refer to points (or loads) rather than surfaces. The general formulation presented in this section applies for both transverse and axial loading.
Displacement and rotation rates are decomposed into elastic and plastic components, that is, The rate of generalized forces is expressed aṡ=̇= where is the force and is the elastic stiffness. The bounding load is constrained through where represents a set of internal variables and̄is the image point. In the one-dimensional case,̄is simply the bounding load magnitude. If we assume that both the current load and the bounding load are centred at the origin, the image point may be expressed through a simple mapping rulē where is the load parameter varying from infinity (when the current load is zero) to unity (when the current load is equal to the image point) and is the ratio between current load and the image point, varying from zero (when the current load is zero) to unity (when the current load is equal to the image point). The relationship between and̄expressed in Equation (4) is referred to as radial mapping. Since the evolution of internal variables is controlled by a function that only depends on the distance between the current load and the image point, Equation (3) may be written as The evolution of plastic displacements is given by the plastic flow rulė where the plastic multiplier is zero during elastic response and greater than zero during plastic response. The hardening rule may be expressed aṡ=̇ (   7) where is the hardening parameter. The consistency conditioṅ̇= and the plastic modulus is expressed as In conventional rate-independent plasticity, the hardening parameter must be defined explicitly in order to obtain the plastic modulus. Here, it is only required that the plastic modulus is (1) a function of (or ) , (2) is infinite during the initial elastic response, (3) vanishes on the bounding point and (4) varies monotonically between the extremes. Several researches 31,32,40 have used the expression to describe monotonic loading for both shallow and deep foundations, where 0 is a calibration parameter. In order to completely describe the cyclic behaviour of foundations, it is necessary to distinguish between virgin loading, reloading and unloading. Chatzigogos et al. 31,32 used the expression to represent a slightly less plastic behaviour during reloading for shallow foundations. Here, is the minimum load parameter obtained during the previous loading steps and is a calibration parameter. However, unloading was described as purely elastic (which is not a realistic assumption for cyclic pile response). Correia 39 and Correia and Pecker 40 simulated the cyclic behaviour of a laterally loaded pile using for both unloading and reloading. Here, 1 is the ratio between the first loading surface and the bounding surface, 2 is the ratio between second loading surface and the virgin loading surface and is a calibration parameter. Equation 14 implies that when current load reaches the virgin loading surface ( 2 = 1), the expression reduces to Equation 12 and the plastic modulus obtained at the last virgin load state is retrieved. The advantage of this formulation is that (1) unloading is inelastic, (2) there is a smooth transition between unloading and reloading and (3) overshooting is avoided at the virgin loading surface. However, this approach produces unrealistically soft behaviour upon partial reloading, especially as the second loading surface expands near the virgin loading surface. In this paper, unloading and reloading are formulated using a different approach, and the formulations depend on the force-displacement relationship. A complete description of the plastic modulus evolution is given in the following sections.

Solution scheme
The numerical scheme is expressed in terms of normalized forces, displacements and stiffness values, i.e., where is the transverse load, is the axial load, is the transverse displacement, is the axial displacement, and are the ultimate (bounding) loads, is the initial, elastic stiffness, indicates a normalized value and is the pile diameter. The algorithm assumes that the initial step in the creation of any loading (virgin loading, unloading or reloading) is purely elastic and that the subsequent steps within the respective loading always contain plastic deformations. To determine the corresponding load and the elastic and plastic displacements, a return mapping algorithm is used. The cutting plane algorithm [60][61][62] has been adopted by researchers 31,32,39,40 in the formulation of macro-elements for both shallow and deep foundations. The same approach will be used here. The main concept of the cutting-plane algorithm is to explicitly integrate the state variables and iterate the solution until a constraining condition is satisfied. Here, this implies using Equations 1, 2, 6, 7 and 10 to obtain the elastic displacement, plastic displacement, load and load multiplier using values from the previous iteration or the last converged step. The consistency condition in Equation 9 is linearised and solved for the plastic multiplier and the state variables are updated. Convergence is achieved when both the displacement residual and the bounding point equation are within a prescribed tolerance value.
The return mapping scheme for the transverse tangent stiffness is shown in Algorithm 1. Here, is the elastic displacement, is the plastic displacement, is the load step number, is the iteration number, subscript represents normalized values, is the transverse load, is the plastic multiplier, is the trial load parameter, is the minimum load parameter obtained during the previous loading steps, is the transverse plastic modulus, is the displacement residual, is the bounding point equation, is the displacement residual tolerance and is the bounding point equation tolerance. The axial tangent stiffness is obtained though a similar procedure.

3.2.1
Plastic modulus evolution for transverse response If , is infinite, , , reduces to the normalized, elastic stiffness , . If , is zero, , , is also zero. In essence, the formulation of a bounding plasticity model condenses into determining the evolution of the plastic modulus between the initial, elastic load (where the tangent stiffness is equal to the elastic stiffness) and the ultimate load (where the tangent F I G U R E 2 Schematic sketch of transverse loading points for reloading and unloading. stiffness is zero). In order to correctly represent the evolution of the plastic modulus within the framework of the presented macro-element, it is necessary to distinguish between virgin loading, unloading and reloading.
Virgin loading is described by the mapping rule such that The plastic modulus is expressed as where 0, , is a calibration parameter. Reloading is described by a different mapping rule. In addition to the current loading point and the image point, it is also necessary to keep track of the virgin loading point and the loading point associated with the last initial unloading and reloading state. Figure 2 shows a schematic representation of loading points for reloading and unloading. Note that we need to separate between reloading in the direction of the last initial unloading point and reloading in the opposite direction of the last initial unloading point. This may be achieved through the parameter where , is the last initial unloading point. To the authors knowledge, this parameter was first introduced by Correia. 39 The image point and the virgin loading point are given by where is the load parameter associated with the virgin loading point. The reloading point is expressed as were , and are the load and load parameter associated with the last initial reloading point. If is positive, reloading is in the same direction as the last initial unloading point. Otherwise, it is in the opposite direction. Note that if is negative, the last initial reloading point and the last initial unloading point represent the same loading point. If 2 is the ratio between the second loading point and the virgin loading point (as defined by Correia 39 and Correia and Pecker 40 ), As mentioned, this formulation may yield somewhat soft response for partial reloading. In this paper, we propose to redefine 2 as the similarity ratio between the second loading point and the virgin loading point span, i.e, F I G U R E 3 Implicit transverse-rotational modification factors.
,1 = 0.5 and ,2 = −0.5 since is equal to when is negative, implying that the plastic modulus is continuous at the transition between unloading and reloading.
Transverse-rotational coupling is implicitly introduced using modification factors for and , , that is, Here, Δ is the current displacement increment, Δ is the current rotation increment, is the displacement at failure for a fixed-head pile, is rotation of a free-head pile at , is the transverse load of a free-head pile at and , is the failure load for pure rotation. Figure 3 shows a plot of the modification factors plotted against the ratio ∕ for typical load ratios. The following should be noted: (a) When ∕ is small, then = = 1. The plastic modulus equals fixed-head conditions. (b) When ∕ = 1, then = = ,1 . The plastic modulus equals free-head conditions. (c) When ∕ is large, then is a negative number restricted by ,2 and is an arbitrarily large, negative number. The plastic modulus equals pure rotation conditions. 2. For −1 ≤ ∕ < 0, then = = 1. The plastic modulus equals fixed-head conditions. It is assumed that when ∕ = −1, transverse loads are mainly caused by transverse displacements. 3. For ∕ < −1: (a) When | ∕ | ≈ 1, then = = 1. The plastic modulus equals fixed-head conditions. (b) When | ∕ | is large, then = ,2 and is a large number. The plastic modulus equals pure rotation conditions.
It is emphasized that the reason for using implicit transverse-rotational coupling is to allow for unrestricted calibration procedures without the need for pre-defined failure surfaces.

3.2.2
Plastic modulus evolution for axial response Similar to transverse loading, it is necessary to distinguish between virgin loading, unloading and reloading. In addition, axially loaded concrete piles behave differently in tension and compression. Virgin loading in tension is described similar to transverse response. In compression, the response is mainly elastic until a cut-off load , is reached. For this range, it is proposed that the plastic modulus varies as , , ( 1 ) = 0, , where is the similarity ratio between the cut-off load and the bounding load, is the similarity ratio between an initial, infinitesimal load and the bounding load and 0, , is a calibration parameter. Once the cut-off load is exceeded, the plastic modulus varies as , , ( 1 ) = 0, , ln Reloading is also defined separately for compression and tension. The similarity ratio 2 is defined as for transverse response, but the variation of the plastic modulus differs in tension and compression. In tension, the plastic modulus varies with Equation (25) (with axial parameters). In compression, it is suggested that the plastic modulus is infinite when reloading is in the direction of the last initial unloading point ( = 1) and , , ( 2 ) = 0, , when reloading in the opposite direction of the last initial unloading point ( = −1). Here, is similarity ratio for the last unloading point and , is a calibration parameter.
Unloading is also defined separately for compression and tension. The similarity ratio 2 is defined as for transverse response. In tension, the plastic modulus varies with Equation (25) (with axial parameters). In compression, unloading is mainly elastic and the plastic modulus is infinitely high. The following features of the formulation should be noted: 1. When virgin loading starts, 1 is low and the plastic modulus expressed in Equation (33) reduces to , , = 0, , ln The initial parameters and are infinite and zero only in a theoretical sense. Numerically, is an arbitrary large number and is an arbitrary small number that need to be predefined in the numerical scheme. Simulations using programming language Python 3 63 show that is a good approximation with respect to desired behaviour and numerical stability. Hence, Equation (36) represents a relatively high plastic modulus, which through Equation (16) (with axial parameters) gives a tangent stiffness value approximately equal to the initial, elastic stiffness. As the load reaches the cut-off load, that is, when 1 is equal to , Equation (33) reduces to Equation (34), implying that the plastic modulus is continuous throughout virgin loading. 2. When unloading in compression, the plastic modulus is given by Equation (36) and the tangent stiffness approximately equals the elastic stiffness. 3. When reloading in compression after unloading from tension ( = −1), and the last unloading point in tension is equal to the bounding load ( = 1), Equation (35) reduces to Equation (25) (with axial parameters). In that case, the plastic modulus is continuous at the transition between unloading and reloading since unloading from tension is also described by Equation (25). However, when the last unloading point in tension is near zero ( is a small number), Equation (35) reduces to Equation (36) and the reloading path equals the unloading path. The calibration parameter , may be considered as a parameter that controls the transition between unloading in tension and reloading in compression between the above-mentioned extremes. When reloading in compression after unloading from compression ( = 1), the plastic modulus is always defined by Equation (36).

Calibration
As mentioned earlier, the developed macro-element does not require pre-defined failure surfaces or other parameters and may be calibrated using any type of non-linear pile-soil model. Vertical pile models may be used for batter piles assuming realistic batter angles up to approximately 15 degrees. There are in total 14 parameters that must be determined. Nine parameters, namely , , , , , , , and , are directly retrieved from the finite element analysis. They may also be taken from closed-form solutions for simple soil profiles. The remaining five parameters ( 0, , , , 0, , , , , ) are obtained by matching the area enclosed by the load path using the macro-element with the corresponding results from the finite element model. The parameters are summarized in Table 1. The given values represent the pile-soil system evaluated in the following sections. The parameters may be determined through the following steps:

Validation
The single pile macro-elements are validated in this section against the finite element model presented in Section 2 for various loading conditions.

Transverse response
Transverse response is shown in Figure 4. Figure 4A shows the cyclic response when rotation is fixed. Virgin loading, unloading and reloading are well defined in a global sense, but there are minor pinching effects caused by the concrete material model that are not fully captured by the macro-element. Figure 4B shows the cyclic response for a combination of in-phase transverse displacements and rotations corresponding to a free-head pile, that is, when ∕ = 1. The stiffness and bounding load reduction are well captured. Figure 4C shows the cyclic response for transverse out-of-phase displacements and rotations. Compared to the case with fixed rotation shown in Figure 4A, the transverse response is only slightly altered. Hence, the assumption that transverse loads are mainly caused by transverse displacements when ∕ = −1 seems reasonable. Figure 4D shows that small displacements with fixed rotation are well captured. Figure 4E shows in-phase small displacements and relatively large rotations. It is observed that transverse forces are negative when transverse displacements are positive. In this case, ∕ is relatively large and rotations also contribute to the generation of transverse forces. When displacements and rotations are in phase, rotations cause negative transverse forces. Figure 4F shows out-of-phase small displacements and rotations. In this case, rotations cause positive transverse forces. The results in Figures 4A-4F show that the implicit transverse-rotational coupling factors expressed in Equations (31) and (32) are able to modify stiffness and bounding load quite accurately for most combinations of displacements and rotations. Figures 4A-4F also show that the transition from reloading to virgin loading is smooth for closed-loop cycles and that the transition from unloading to reloading (point of zero loading) is smooth. Figures 4D and 4E demonstrate that the macro-element is capable of capturing cyclic behaviour for small transverse displacements. However, it is emphasized that the macro-element is not able to simulate frequency-dependent radiation damping, which becomes more important with decreasing displacement amplitudes. Figure 4G shows one partial loading cycle. In addition to comparison with the finite element model, the macro-element is demonstrated using 2 defined by Equations (23) and (28). In this case, it is clear that 2 defined by Equations (24) and (29) produce more accurate behaviour, especially as the load approaches the virgin loading point. Figure 4G also shows that even though the plastic modulus is discontinuous in transition from reloading to virgin loading, the behaviour is still very close to the finite element model. This behaviour is supported by the findings in The Pile Soil Analysis (PISA) Project, 44,64 where design procedures for OWT monopiles were addressed through numerous medium scale field tests. Figure 4H shows multiple partial loading cycles. The macro-element is generally able to capture this behaviour fairly well, but is should be emphasized that overshooting is only avoided at the virgin loading point. This is shown in the blown-up segment in Figure 4H. However, this error mainly occurs in the vicinity of the virgin loading point with small magnitude. As for the case with one partial loading cycle, using 2 defined by Equations (24) and (29) produces better results compared to Equations (23) and (28.)

Axial response
Axial response is shown in Figure 5. Figures 5A and 5D show that cyclic loading with large amplitudes is well simulated. Note that the macro-element is able to capture complex features such as (1) the highly elastic response for virgin loading and unloading and (2) the abrupt stiffness change passing the cut-off load in compression. Figures 5B and 5E show that the behaviour under small displacements is also well captured. As for transverse response, it should be noted that the macro-element is not able to simulate frequency-dependent radiation damping. Figures 5C and 5F show that partial loading is fairly well represented, admittedly with some inaccuracy. Most importantly, Figure 5F shows that Equation (35) is able to globally capture the transition between unloading in tension and reloading in compression for a variety of unloading points. The observed axial behaviour in compression is supported by experimental tests reported in the literature. Chen et al. 65 performed a number of cyclic, axial loading tests on pre-stressed concrete pipe piles in soft clay. The load-displacement curves for quasi-static compression tests showed similar trends for virgin loading and unloading as observed in this study. Drbe and El Naggar 66 performed full-scale load tests on micro-piles in cohesive soil. Quasi-static tests for loading, unloading and partial loading in compression showed the same tendency as the results obtained in this study.

Stiffness matrix assembly
The single pile formulations for the individual load-displacement relationships has been presented in the previous sections. This section provides a procedure for establishing a three degree-of-freedom stiffness matrix based on the single pile formulations. The macro-element is validated for a variety of pile configurations, batter angles and soil profiles. It is demonstrated how the macro-element may be employed to obtain linear-equivalent properties. With reference to Figure 6, the displacement increment of the single pile is expressed as where is the displacement vector of pile in local coordinates, is the coordinate transformation matrix and is the displacement vector of pile in global coordinates. The tangent stiffness matrix of a single pile in global coordinates is obtained through a coordinate transformation procedure and is expressed as where , is the tangent stiffness matrix of a single pile in local coordinates and , , and , , are the tangent stiffness values returned by the respective single pile formulations. The pile group stiffness matrix is expressed as The matrix entries in are obtained by enforcing unit displacements and rotations (in turn) and solving the equilibrium equations, that is, where is the distance from the global node to the pile-node and is the number of piles. We assume that there is no contact between the pile cap and the soil, that is, all the forces and moments from the pile cap are transferred through the piles. Note that (1) is set negative when the corresponding pile is located on the right hand side of the global node, (2) the pile-cap rotation is resisted by vertical pile-head loads only and (3) the local nodes are internal element nodes that are not a part of the global solution.

Validation
A 2×1 pile group with vertical piles and 0 ∕ = 5 is subjected to harmonic, in-phase horizontal displacements and rotations. The results are shown in Figure 7. Figure 7A shows the results for large displacements and small rotations. Figures 7B and 7C show the same results, but with increasing rotation. It is clear that the implicit transverse-rotational coupling modifies horizontal stiffness and bounding load quite accurately in all three cases. The moment-rotation relationships match the finite element model fairly well, but it is evident that the accuracy decreases for small rotations. The same analyses are carried out for a 2×1 pile group with batter piles ( 0 ∕ = 5 and = 15 ). The results are shown in Figure 8. As for the vertical pile group, the macro-element shows good agreement with the finite element model. The moment-rotation relationships are particularly well captured, but the horizontal load-displacement relationship is somewhat overestimated. Interestingly, the horizontal stiffness of a batter pile group is practically unaffected by rotation. Since the opposite applies for vertical pile groups, the results indicate that batter piles substantially increase horizontal stiffness, especially at large rotations. Next, a 3×3 pile group with both vertical and batter piles ( 0 ∕ = 5 and = 15 ) is subjected to various combinations of horizontal displacements and rotations. The results are shown in Figure 9. Figure 9A shows the results for in-phase large displacements and small rotations. As for the 2×1 pile group with batter piles, the horizontal load-displacement relationship is somewhat overestimated. The same accuracy is observed in Figure 9B for in-phase large displacements and large rotations. Figures 9C, 9D and 9E show that the macro-element is able to capture in-phase small displacements and large rotations, out-of-phase large displacements and small rotations and out-of-phase large displacements and large rotations. Figure 9F shows that partial loading is captured with approximately the same accuracy as full cycle loading.
One of the limiting assumptions governing the macro-element formulation is that pile-soil-pile interaction is neglected. Figure 10A shows the results for a 3×3 pile group with 0 ∕ = 3 in the direction of loading. The finite element model shows only slightly softer horizontal response compared to the pile groups with 0 ∕ = 5 shown in Figure 9B. These results indicate that neglecting pile-soil-pile interaction may be reasonable for pile groups with realistic 0 ∕ -values and soft soil. It is, however, expected that the accuracy will decrease as the soil stiffness increases. In addition, it was suspected that neglecting pile-head moment would decrease the moment-rotation accuracy for close pile spacing. However, rotation is well captured in this case, and the difference between 0 ∕ = 3 and 0 ∕ = 5 is evident in Figures 10A and 9B.
Since the numerical scheme is based on de-coupled, single pile response, the macro-element should not be restricted to symmetric configurations. Figure 10B shows the results for a 3×2 pile group with asymmetric configuration of batter angles in the direction of loading. The macro-element is in good agreement with the finite element model.
All analyses thus far are performed using the (approximately) linearly varying soil profile shown in Figure 1B. In order to validate the macro-element for different soil profiles, the 3×3 pile group is analysed in homogeneous soil corresponding to the middle layer in the linearly varying profile. The results in Figure 10C show that the macro-element matches the finite element quite well. The results presented herein indicate that even though the macro-element may be somewhat inaccurate compared to a rigorous finite element model, it is highly capable of capturing trends associated with pile group configuration, batter angle, pile spacing and soil profile. In practical engineering, where actual response values are approached in a broader sense, it is often more important to evaluate trends rather than exact values. The macro-element is therefore particularly suitable in practical design situations or as an efficient tool in parametric analysis for both practical and academic purposes.

Equivalent linear properties
In addition to non-linear time history analysis, the macro-element may be used to efficiently estimate equivalent linear properties of a pile group. The equivalent linear pile group impedance in complex form is expressed as where is the displacement-dependent secant stiffness and is the equivalent viscous damping ratio. The equivalent viscous damping may be estimated using the area based approach first suggested by Jacobsen, 67 that is, where is the energy dissipated in one cycle during the steady state response and is the peak energy during one cycle. The equivalent linear approach is schematically shown in Figure 11A for the 3×3 pile group subjected to harmonic displacements with fixed rotation in linearly varying soil.
The macro-element allows for fast parametric analysis of secant stiffness and damping ratios for a variety of parameters such as pile group configuration, batter angle and pile spacing. As an example, Figure 11B and 11C show the secant stiffness and damping ratios for a variety of batter angles as a function of displacement amplitude. In this case, it is clear that secant stiffness is highly dependent on batter angle. The damping ratio, however, is rather unaffected by batter angle. The ability to efficiently asses such values is particularly useful in preliminary design or as means of evaluating the effect of change in later design stages. In addition, governing codes and general design rules are often expressed in terms of simplified, linear values, which in most cases are not easily retrieved using rigorous finite element tools.

CONCLUSION
A novel macro-element for vertical and batter pile groups has been presented. The numerical scheme is based on decoupled, single pile response without any requirements for pre-defined failure surfaces or other parameters. Although practical, such simplified formulations are bound to limited validity. First, the accuracy is expected to decrease for small displacement since the formulation neglects radiation damping, which becomes more important for small-strain soil deformations. Second, even though the macro-element performs well within the realistic range of pile spacings in soft soil, further studies on the performance in stiffer soils are needed. Third, axial and transverse response are de-coupled, which inherently introduces an error related to the bounding loads and restricts the macro-element to long piles. Nevertheless, it has been demonstrated that the macro-element is capable of capturing trends associated with pile group configuration, batter angle, pile spacing and soil profile.

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

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.

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.