Numerical model for the dynamic and seismic analysis of pile-supported structures with a meshless integral representation of the layered soil

This paper presents a three–dimensional linear numerical model for the dynamic and seismic analysis of pile-supported structures that allows to represent simultaneously the structures, pile foundations, soil profile and incident seismic waves and that, therefore, takes directly into account structure–pile–soil interaction. The use of advanced Green’s functions to model the dynamic behaviour of layered soils, not only leads to a very compact representation of the problem and a simplification in the preparation of the data files (no meshes are needed for the soil), but also allows to take into account arbitrarily complex soil profiles and problems with large numbers of elements. The seismic excitation is implemented as incident planar body waves (P or S) propagating through the layered soil from an infinitely–distant source and impinging on the site with any generic angle of incidence. The response of the system is evaluated in the frequency domain, and seismic results in time domain are then obtained using the frequency–domain method through the use of the Fast Fourier Transform. An application example using a pile-supported structure is presented in order to illustrate the capabilities of the model. Piles and columns are modelled through Timoshenko beam elements, and slabs, pile caps and shear walls are modelled using shell finite elements, so that the real flexibility of all elements can be rigorously taken into account. This example is also used to explore the influence of soil profile and angle of incidence on different variables of interest in earthquake engineering.

An adequate modelling of the wave propagation phenomena through the soil and of the pile-soil-pile interaction effects is needed when studying the dynamic and seismic behaviours of pile foundations and of pile-supported structures. Moreover, it is also of interest to be able to model in detail the variability of the soil properties with depth and the incoming seismic waves (that are also affected by the soil profile). However, in many cases, taking into account all these elements might imply the need of building complex models, or might even prove to be too costly in computational terms, leading to simplified models. This idea justifies the efforts in the line of developing a rigorous numerical model for this problem allowing the study of configurations with large number of foundation elements and complex soil profiles with reasonable low computational requirements Among the different numerical methods that can be used for the analysis of the interaction between piles and soil, the boundary element (BE) method is one that stands out for its capability to treat the radiation of energy through the unbounded media implicitly in its formulation. However, the significant increase in terms of computational effort and memory usage needed to model a large number of layers, or the difficulty to represent continuously varying soil properties, are significant disadvantages that limit the applicability of the boundary element method to some problems, unless specific advanced fundamental solutions are employed to model layered soils without the need of discretizing the layer interfaces.
For this reason, Álamo et al. (2016) proposed a numerical model for pile foundations in which, by employing advanced Green's functions for modelling the dynamic response of the layered unbounded viscoelastic soil medium together with a Timoshenko beam formulation for modelling the dynamic response of piles, the integral equations that describe the dynamic behaviour of the soil-foundation system are significantly simplified and produce a very compact numerical approach for the problem. No meshes are needed to represent the geometry of the soil or its profile, and only few degrees of freedom are required to represent the piles, so that problems with complex soil profiles and large number of pile foundations can be studied in detail and including all media, with a limited number of degrees of freedom.
In the work at hand, this soil-pile integral model previously presented by Álamo et al. (2016) is further developed by coupling the pile foundations with finite element structures, as well as including seismic and superficial excitations. Thus, the proposed model can be used both for a direct approach of the soil-pile-structure problem, and for studying the behaviour of pile foundations in order to use this information in substructuring procedures. The hypotheses and formulation of each part of the model are detailed in Sect. 2. Then, results are compared to the ones of a previous BE-FE model in Sect. 5.2, followed by Sect. 5.3 in which an example problem is analysed in order to highlight the capabilities of the proposed tool in earthquake engineering. Finally, Sect. 6 summarizes the main conclusions of the work.

Model description
The model proposed herein (see Fig. 1) is a three-dimensional, linear numerical model formulated in the frequency-domain for the study of the response of pile foundations and pile-supported structures through a direct approach. The system can be subjected to planar wave-fronts (type P or S) propagating with a generic direction within the soil. The seismic response in time domain under a specific earthquake is obtained using the frequency-domain method (Chopra 2017) through the use of the Fast Fourier Transform. Other types of loads, such as prescribed displacements or forces acting at different points of the structures/piles, or loads acting on the free-surface of the soil domain, can also be considered. Material damping of hysteretic type is adopted through the definition of hysteretic damping coefficients ( ) and the corresponding complex properties (see e.g. Christensen (1982)).
First, the soil dynamic behaviour (Sect. 2.1) is represented by a formulation based on the integral reciprocity theorem, the use of advanced fundamental solutions and the treatment of piles as load lines acting within the soil. Then, the way in which the piles are introduced in the model through beam finite elements is discussed in Sect. 2.2. Superstructures (Sect. 2.3) are modelled through FE with the combination of beam and plates elements. The coupling between the superstructure and pile elements is made by imposing compatibility (in terms of nodal displacements and rotations) and equilibrium (through pile-structure coupling forces) conditions, and the solution of the whole coupled system can be found from the resolution of a linear system of equations involving only unknowns related to piles and structures (Sect. 2.4).

Soil equations
The reciprocity theorem in elastodynamics (see Wheeler and Sternberg (1968)) relates two different elastodynamic states denoted as S(u j , ij , b j ; , Ω) and S * (u * j , * ij , b * j ; , Ω) , that satisfy the Navier equations in the domain Ω . Assuming harmonic and zero-initial conditions, the integral equation of this theorem can be expressed through Einstein's notation as: where u j and u * j are the displacements at any point of the domain Ω ; p j and p * j are the tractions acting over the boundary Γ = Ω that are compatible with the stress tensors ij and * ij , respectively; and b j and b * j are the body forces acting inside the domain Ω , all for each of the two states mentioned above.
The unknown state S corresponds to the actual problem to be solved, while the state S * is a known state corresponding to a point source solution that is usually referred to as fundamental solution and whose proper definition can significantly simplify the formulation of the problem. All variables related to this known state are denoted by the superindex*. For the proposed model, the Green's functions developed by Pak and Guzina (2002) are used as fundamental solution. These Green's functions represent the response at any point of a layered half space to a test harmonic unit point load applied at a certain collocation point (let's denote it as point : ) in each one of the three directions of the space. For this reason, the known displacements, tractions and body forces corresponding to the known state S * in Eq. (1) must be now understood as tensors u * ij , p * ij and b * ij , respectively, where index i denotes the direction of the unit test load, and where p * ij = 0 at Γ , and b * ij = ( ) ij , being ( ) the Dirac's delta and ij is the Kronecker's delta. Thus, the domain integral of the left-hand side of Eq. (1) can be computed as the displacement at such collocation point u i . These Green's functions are obtained by solving the Navier's equations in terms of displacement potentials (see Pak (1987)) in a transformed cylindrical coordinate system that uses Fourier and Hankel transforms for the angular and radial coordinates, respectively. The formulation includes particular propagation matrices to relate the response of the different layers. These matrices do not contain any unbounded exponential terms in order to avoid numerical instabilities. A particular integration procedure, as described by Pak and Guzina (2002) and Pak (1999, 2001) is employed to evaluate the inverse Hankel transform, which has been modified for low-frequencies based on the single layer problem (Martínez-Castro and Gallego (2007)). Note that, as these Green's functions are obtained assuming a layered half space domain, the boundary conditions of the free-surface and layer interfaces are intrinsically satisfied by the fundamental solution. Therefore, the contour integral of the left-hand side of Eq. (1) vanishes as the fundamental solution satisfies the free-surface conditions (i.e., p * j = 0 at Γ). From the point of view of the soil equations, the effect of the piles is modelled through load lines Γ p acting inside the soil domain. The interaction forces between soil and pile are reduced to a distributed force q j along a line located at the axes of the piles and representing the resultant of the tractions acting over the soil at the soil-pile interface (pile shaft) as a function of depth and direction j.
In the absence of external forces acting over the free-surface in the unknown state S , the contour integral of the right-hand side of Eq. (1) would also vanish ( p j = 0 at Γ ). On the other hand, the body forces of this unknown state are limited to the distributed tractions q i acting over the soil due to the soil-pile interaction along the soil-pile interfaces Γ p , which allows to simplify the domain integral of the right-hand side of Eq.
(1) to a line integral along Γ p . Thus, the integral equation that could be used to model the soil response in the absence of incident seismic waves and external actions over the free surfaces would be the following compact expression: which involves only the pile-soil interaction tractions, the displacements along the piles axes, the fundamental solution in terms of displacements and line integrals to be evaluated only along the piles axes. Due to the characteristics of the fundamental solution, this integral becomes singular when the collocation point belongs to the integration element. In order to avoid this singularity, a non-nodal collocation strategy as described in Álamo et al. (2016) and Álamo (2018) is adopted.
On the other hand, the seismic excitation can be considered through the classic decomposition of the total field into the incident and scattered fields and rewriting the integral equation of the soil in terms of the variables of the scattered field (see e.g. Domínguez (1993) and Padrón et al. (2008)). The travelling wavefronts determine the incident field of displacements u I j and stresses I ij at each point of the layered half space domain. The incident field can be analytically obtained through the resolution of wave propagation problems corresponding to different types of waves (Section 4 briefly presents their expressions for the wave types that can be considered in the proposed model). In this case, in which the soil profile is already taken into account when computing the incident field, the scattered fields arise from the presence of the piles in the soil. The stiffness of the piles makes them oppose the motion of the soil, creating this additional field that propagates from the piles to the surrounding soil. The dynamic response of pile foundations and superstructure will also generate part of this field. Therefore, when the system is excited by an incident field, Eq.
(2) and its variables have only meaning for the field scattered by the piles (and the superstructures), and it should be rewritten in terms of the total field as: where u i here represents the (unknown) total displacement constituted by the superposition of the incident and scattered fields at the collocation point; while u I i is the (known) incident field displacement at the same point, formulated in the absence of piles and structure and incorporating the boundary conditions of free surface and layer interfaces of soil (see Sect. 4).
After pile discretization, and by applying Eq.
(3) at every pile node, the system of equations corresponding to the soil can be written as: where ̄ and ̄ are the vectors of nodal displacements and soil-pile interaction tractions, ̄ is the vector containing the incident field displacements evaluated at the pile nodes, and is the influence matrix, that encapsulates the dynamic response of the soil and includes its material and radiation damping. This influence matrix is obtained by numerical integration of the fundamental solution times the shape functions of the interaction tractions. Details about the discretization (elements and shape functions) and numerical evaluation of the integrals are given in Sect. 3.
Once the system response is known, in a post-processing stage Eq. (3) can be applied at any internal point of the soil domain in order to compute its displacement vector as: where the superindex in the influence matrix indicates that the collocation point is the point , while are the incident field displacements evaluated at this point of the soil.

Pile equations
Beam finite elements are used to model pile lateral and axial behaviour. On the contrary, the torsional mode of piles is neglected as the soil-pile model does not contemplate any interaction mechanism for it. Following the typical assembly procedure, the finite element system of equations corresponding to the piles can be written as: where ̄ is the vector of pile nodal displacements and rotations, ̃ = − 2 contains the global stiffness (including hysteretic damping) = Re[ ](1 + 2i ) and mass matrices of the piles, and the right-hand side of the equation contains the forces acting over the piles. These forces are divided into two components. The vector ̄ contains the forces (and moments) acting at the pile head, which are assembled into the correct equations of the global system through index matrix . The distributed soil-pile interaction tractions acting along the pile shaft are defined trough their nodal values ̄ and converted to nodal loads through the global matrix (see Álamo 2018 for details). In the general scenario, where the pile head forces are produced by the coupling with the superstructure, all of the aforementioned vector variables are unknowns of the problem. Thus, the pile system of equations should be written as: Any pile caps would be modelled as structural shell elements, as described below. This way, the flexibility of these elements of the foundation can also be taken into account.

Structure equations
The analysis of the superstructures is conducted by a FE representation with unidimensional elements (beams) and bidimensional elements (shells) that can be freely combined. The only restriction of the structure model is that the different elements are connected through their nodes. In addition to these elements, the existence of punctual masses or moments of inertia on the structural nodes is also allowed through their inclusion into the corresponding terms of the structural mass matrix.
By assembling the elemental stiffness and mass matrices of all of the structural components, the FE system of equations corresponding to the structures can be written as: where ̄ is the vector of structural nodal displacements and rotations, ̃ = − 2 contains the global stiffness (including hysteretic damping) = Re[ ](1 + 2i s ) and mass matrices of the structural elements, ̄ is the vector of prescribed nodal forces and moments acting over the structures, and ̄ is the vector of forces and moments acting over the structures due to the coupling with the pile heads. These coupling forces are correctly introduced into the corresponding equilibrium equations with the aid of indexing matrix .
Rearranging Eq. (8) in order to separate the known and unknown variables, the system of equations corresponding to the structural elements results in:

Coupling equations
The coupling between soil and pile variables is made by imposing compatibility and equilibrium conditions in terms of displacements and interaction tractions, respectively, at the soil-pile interface. This assumption implies that no sliding or gaping phenomena can be directly considered in the analyses, although their effects can be approximated through equivalent linear models (González et al. 2020). However, the considered fullcontact assumption is able to reproduce the foundation response under small amplitude vibrations and also the physical mechanisms of the interplay between piles and soil freefield response. For easing the coupling, the same discretization is used for the pile elements and soil load lines. Thus, soil-pile coupling equations can be simply written as: where is an auxiliary matrix that extracts the displacement terms of the pile displacements and rotation vector and that is also needed for non-nodal collocation strategies (see Álamo (2018)).
On the other side, piles and structures are assumed to be coupled at the pile heads. In the proposed model, the coupling is directly made between the pile head nodes and the structural nodes located at the same point of space. Compatibility between the corresponding nodal displacements and rotations, and equilibrium between the pile head forces and structure-foundation forces can be expressed trough: where and are the matrices defining the compatibility equations between the corresponding nodes of structures and piles. No direct contact has been considered between the structure and the soil, i.e, the interaction with the soil can only be modelled through the piles, and any structural element defined over the free-surface is assumed to be not in contact with it.
Introducing the coupling equations from (10) to (13) into the governing equations of soil (4), piles (7), and structures (9); the global system of equations can be written as: which is a system of equation written in terms of variables of pile and structure only, without specific unknowns of the soil.

Loads on free-surface
Apart from the seismic loads considered above, prescribed distributed harmonic forces can be imposed over a portion of the free-surface. Let p ext j be the load acting, in the direction j, over a part Γ ext of the soil boundary Γ , i.e. free-surface. For generality's sake, p ext j can represent different load types, such as a point, line or surface load. The contour integral of the left-side of Eq. (1) does not longer vanish but can be reduced to Γ ext : Thus, the integral equation of the soil with the inclusion of the external load at the freesurface results in: where two-and four-noded elements (see Fig. 2) are used for discretizing the line and surface loads, respectively, in order to evaluate the numerical integral which, in the case of a point load, simplifies to a single multiplication.
Considering the pile and the external load discretizations and applying Eq. (16) to all pile nodes, the soil system of equations including the loads on the free-surface is obtained as: where ̄ ext is the vector containing the nodal values that defines the external load on the free-surface and ext is the influence matrix obtained from the numerical integration of the fundamental solution times the proper shape functions used to define the external load inside each load element. Note that the vector ext̄ ext is known and forms part of the righthand side once the equations are rearranged into the final system of Eq. (14).

Speed-up procedure for the computation of influence matrices
The accurate computation of the terms of the soil influence matrix requires a large number of evaluations of the fundamental solution. This process is computationally very time consuming due to the iterative and complicated nature of the used Green's functions for the layered half space. Of course, this issue becomes more important as the number of piles or their discretization increase and, also, for high frequencies or complex soil profiles.
However, as the value of the displacements of the Green's functions only depends on the relative distance between the collocation and observation point and on their depths with respect to the free-surface, some terms of the influence matrix are duplicated or, at least, equivalent through an in-plane rotation. As the nodes of each pile are located at different depths, the equivalence between different sub-matrices is established at a pile scale. Let The reduction in computational effort and time can be very significant in the case of regular pile groups, and increases with the number of piles in the group. In square regular groups of N × N piles, the reduction in computational time obtained with this strategy has been found to be above 40%, 70%, 80%, 90% and 95% for N = 2, 3, 4, 6 and 10, respectively.

Element types
The unidimensional elements are modelled by using 2-noded beam elements. Structural elements (e.g., columns) present 6 degrees of freedom per node (3 displacements and 3 rotations), while pile elements have 5 degrees of freedom per node (3 displacements and 2 rotations) as their torsional mode is neglected. Cubic and quadratic shape functions, as described by Friedman and Kosmatka (1993), that satisfy the static governing equation of the Timoshenko's beam, are used in order to define the lateral displacements and rotations of the element, respectively. Linear shape functions are used to represent the axial variables.
On the other hand, building slabs and pile caps are modelled through the Mixed Interpolation of Tensorial Components (MITC) shell elements (see Bucalem and Bathe (1993); Dvorkin and Bathe (1984)). These elements, based on the degenerated solid approach (see Ahmad et al. (1970)), are chosen due to its predictive capacity and robustness, as their formulation avoids the typical issues of locking or spurious modes (see Yang et al. (2000)) of shell elements. In particular, 4-noded and 9-noded quadrilateral elements are considered. Each node has 6 degrees of freedom corresponding to the three displacements and three rotations in space. The choice of the 6 (three global rotations) degrees of freedom over the 5 (two local rotations) ones is made in order to ease the coupling between the different structural elements.

Expressions for the computation of the incident fields
The seismic excitation is modelled as a wavefront of body waves originated by an infinitely-distant source. This wavefront is assumed to propagate through the half space layer inside the x − z plane with a known (yet arbitrary) angle of propagation. Once this incident wavefront reaches the first layer interface, reflection and transmission phenomena take place, producing different waves that either propagates upwards to the free-surface or are reflected back to the unbounded media. Figure 4 shows a sketch of the problem, together with the layer numbering and coordinate system of reference that is used along this section in order to obtain the expressions of the incident fields. Here, h n is the layer n thickness; c n s and c n p are their Shear and Primary wave propagation velocities; n is its shear elastic constant; and n SH , n SV and n P are the angles of propagation of each wave type. In the following, only the main aspects of the procedure are presented, for a step-by-step explanation see (Wolf 1985).
The trigonometric functions of each angle of propagation n W corresponding to the waves of type W in the layer n will be denoted as: −1 the imaginary unity. This definition of t n W results in the correct sign for the case of SV waves with an angle of propagation below the critical one (see e.g. Wolf (1985)).
For each layer n, the displacements and stresses produced by the incident field can be obtained as: where z is the local vertical coordinate of the layer (see Fig. 4) and k = ∕c is the wavenumber obtained from the corresponding apparent velocity: The apparent velocity indicates the speed at which the wave front travels along the horizontal direction and its value is the same for all of the soil layers in order to satisfy the compatibility condition in the x direction for each soil interface. Thus, from Eq. () it is For the out-of-plane motion problem, the depth-varying amplitudes for the displacements and stresses at each layer can be obtained from the amplitudes of the incident and reflected SH waves following: where For the in-plane problem, the expressions between the non-zero displacements and stresses and the amplitudes of the upwards and downwards SV and P waves in each layer are: where: The frequency-dependent amplitudes of all waves involved in the problem can be computed after solving the linear system of equations obtained by imposing the following boundary conditions: -free-surface conditions: -compatibility conditions at each interface n between layers n and n + 1 : -incident wave type at the deepest layer N + 1 (half space): SH waves (out-of-plane motion): Once the wave amplitudes are known, the free-field response of the soil (incident field) can be computed by using either Eq. (23) or Eq. (25) depending on the wave type.

Results
In order to illustrate the capabilities of the model described above, this section presents results corresponding to the seismic response of a pile-supported building founded on a stratified soil and subjected to different earthquakes generated by trains of S-waves that impinge the site vertically or with a certain angle of incidence . Section 5.3 analyzes a problem that exploits many of the capatibilities of this new model through the analysis of a pile-supported structure, subjected to inclined seismic SV waves, and founded in a soil deposit with properties that change continuously with depth and that overlays an elastic bedrock. However, before that, the code is validated in Sect. 5.2 through comparison against results of a previously published BE-FE model by Padrón et al. (2011) for a simpler base configuration. Firstly, the problem studied in both sections is described below. Figure 5 presents a sketch of the problem at hand of a pile-supported building subject to seismic SV waves. The structure is founded on a soil deposit overlying a stiff homogeneous half-space. Two different profiles will be considered for the soil deposit: a homogeneous profile (type A) with a constant Young's modulus E A = 0.28 GPa ; and a continuously inhomogeneous profile (type B) with a Young's modulus that increases linearly with depth, from E B = 0.1E A at ground surface to E B = 1.9E A at the bottom of the deposit. Both deposits present the same mean value of the Young's modulus. All relevant properties are described in Fig. 5. The properties of the five-storey RC building and its square 4 × 4 pile foundation are presented in Fig. 6. It is worth noting that the BE-FE code considers rigid slabs, with inertial properties lumped at the slab geometrical centres, while the proposed Similarly, the former model assumes massless columns, while the second model assumes a density for the concrete in the columns coincident with that in the piles.

Verification results
The implementation of the proposed model is validated through comparison of relevant variables against those computed using a previously presented BE-FE code by Padrón et al. (2011). Soil deposit type A and vertically incident S waves are considered herein. For the sake of brevity, only inter-storey drifts for floors 1, 3 and 5, and corner pile head bending moments, shear forces and axial forces are compared, although the rest of variables show also an excellent agreement. Figure 7 presents the comparisons in the frequency domain, with the aim of analysing not only the level of agreement but also of being able to detect any frequency-dependent discrepancies. Inter-storey drifts are plotted in the left plots, and pile head forces (bending moment M, shear force V and axial force N) are represented on the right. All functions are plotted in the range 0-10 Hz with respect to the free field ground motion ( u ff ). The proposed model agrees very closely with the previous BE-FE code taken as reference here if the same simplifying assumptions are considered in both models (rigid slabs and lumped inertia properties). It is worth remembering here that such BE-FE code requires a much larger number of degrees of freedom to analyse a problem like this due to the need of discretizing the ground surface and any soil interfaces. For this reason, computing times and memory requirements are reduced when using the proposed model, and the analysis of problems with large numbers of interfaces becomes prohibitive with the former BE-FE code which, at the same Fig. 6 Two-dimensional sketch of considered pile-supported structures time, is not able to represent layers with continuously varying elastic properties. All these issues were the incentive for the development of the present model.
On the other hand, the proposed model allows to study how the dynamic response of the system changes when the slab flexibility and distributed inertia properties are assumed. As shown in Fig. 7, the shell behaviour of the structural slabs significantly reduces the stiffness of the system, shifting the peaks corresponding to its natural frequencies to lower values, and, in general, reducing the amplitude of the inter-storey drifts. The shell behaviour of the slabs and the variations in the soil-structure system natural frequencies are illustrated by Fig. 8, where the structural deformation is presented at the first three natural frequencies. A colour scale is used in order to highlight the distribution of the vertical displacements through the slabs. Regarding to the effects on the corner pile head forces, the axial and shear forces are reduced, while its bending moment is, at small frequencies, amplified if the slab shell behaviour is considered. These phenomena are strongly influenced by the significant flexibility assumed for the pile cap in this illustrative example. If a rigid cap is considered instead, the pile head bending moment becomes clearly proportional to the pile head shear force (as explained  (2021)), with a significant decrease in the magnitude of the head bending moments and an increase in the magnitude of the head shear forces.

Application results
Finally, the results obtained for the problem described in Sect. 5.1 considering both types of soil deposit and inclined SV incindent seismic waves are presented here assuming that the site is impinged by five different earthquakes. The accelerograms are extracted from   (2021), and correspond to seismic events measured in stations located over soils with mean shear wave velocities V s,30 close to 240 m/s. For reproducibility's sake, Table 1 identifies the used accelerograms, indicating the Record Sequence Number (RSN) of the database and the horizontal component used for the analyses, as well as information about the earthquake event and measuring station. These acceleration signals are assumed to correspond to the free-field motion at surface level, independently of the angle of incidence or soil profile. This assumption is made in order to illustrate how much the uncertainties about soil and wave propagation mechanisms can alter the structural response even if the free-field motion at surface level is known. In order to compare the results of all signals, they have been scaled in order to present the same maximum ground acceleration (0.2 g). Figure 9 presents the normalized acceleration spectra of the selected accelerograms. Figure 10 presents the envelopes of seismic bending moments along the piles for three different angles of incident of the SV waves: = 50 , 75 and 90 o , corresponding = 90 o to vertical incidence. Results are presented for selected piles in the group, as illustrated in the plots. The direction of propagation of the waves is also depicted in the sketches. Results corresponding to soil deposit types A and B are plotted in black and red curves, respectively, with dashed lines used for the envelopes obtained for each individual earthquake, and solid lines for the average envelope. The expected local peak bending at the interface (see e.g. Dezi et al. (2010); Sica et al. (2011);) is observed, with a more significant value in the case of the homogeneous soil deposit (case in which the stiffness contrast at the interface is largest). The dispersion among different accelerograms is small, with the exception of the response around the pile head, especially for the corner piles and = 50 o . These corner piles present a significantly higher bending moment than the central ones, denoting a higher contribution of the inertial-induced moments.
In order to analyse this effect, Fig. 11 compares the obtained average maximum bending moments along the piles (solid lines) with the results obtained for the kinematic only problem (dashed lines), i.e. pile foundation without the structure and assuming a massless flexible cap. The comparison clearly shows that the pile bending moment is determined by the kinematic problem once a certain depth is reached. Near the pile head, the combination of inertial and kinematic loads (being not necessarily in phase) modify the kinematic response of the foundation, being this effect more important at the corner piles. The bending profiles of piles affected by the oscillation of the superstructure Fig. 9 Response spectra of the used seismic signals are larger for profile B (variable properties) than for profile A (homogeneous). This fact qualitatively agrees with the pile active lengths corresponding to the studied profiles (Karatzia and Mylonakis 2016): 5.4 m (profile A) and 6.1 m (profile B). On the other hand, the influence of the angle of incidence in the maximum kinematic bending moments is not significant, with the exception of the corner pile at surface level. For these piles, as the angle of incidence gets lower, higher bending is produced at several meters below the surface. Comparing the kinematic-only results for the two studied pile positions, a noticeable group effect near the pile head is found, which increases as the angle of incidence of the waves becomes smaller. It is worth highlighting that this group effect is mainly due to the presence of the flexible massless pile cap that is considered in the kinematic only case.
The influence of the angle of incidence of the SV waves on the response of the soilfoundation-structural system is closely studied in Fig. 12, which presents the evolution of the average maximum seismic pile head bending moments with the angle of incidence for all piles in the group in the range = 50 o to 90 o . The greatest variations are found for the rear central and external piles in a homogeneous deposit (solid blue curves in top left and bottom left plots, respectively) with the maximum seismic head bending moments increasing by 65% from = 90 o to = 50 o . In the rear corner piles (right bottom plot), where the absolute maximum pile bending moments are found, the bending   decreases by around 7% in that case. In the foundations embedded in the inhomogeneous soil deposit, bending moments tend to be larger (except for the corner piles) but the influence of the angle of incidence is smaller (probably because in this case, the waves tend to the vertical direction due to the variation in the soil properties).
On the other hand, maximum inter-storey drifts are presented in Fig. 13 for = 50 , 75 and 90 o . The height of the five slabs at z = 4 , 8, 12, 16 and 20 m is indicated in the plots. As usual, the maximum drift is found for the first floor. The influence of the angle of incidence on this variable is generally not large, with average differences with below 10%, although the maximum difference can be found in the inter-storey drift at the first floor, with an increase of 40% from = 50 o (19.44 mm) to = 75 o (27.26 mm) in the case of the average response for the soil deposit type B (and with similar values for type A).

Conclusions
A new model for the dynamic and seismic analysis of pile-supported buildings including structure-foundation-soil interaction, seismic incident waves and the possibility of modelling complex soil profiles without an excessive increase in computational effort has been presented in this paper. Thus, this model allows the direct analysis of the coupled problem (including structures, foundations, soil and seismic action) without the need of substructuring methodologies and with no further simplifications of the problem beyond the implicit simplifying assumptions of linear-elastic behaviour of soils and structural elements, and bonded contact conditions between piles and soil. The model is formulated in the frequency domain, and seismic response in time domain for a certain accelerogram is then obtained using the frequency-domain method through the use of the Fast Fourier Transform. In order to maximise the flexibility of the code, piles and columns are modelled using Timoshenko beam finite elements, while structural slabs, pile caps and shear walls are modelled using shell finite elements. The more significant advantages or the model are its flexibility, the possibility of modelling complex problems with large numbers of elements, and the low effort required to prepare the data files, as no meshes are needed to define the soil surface or interfaces. Furthermore, the use of a direct approach of the problem instead of, e.g., a substructuring methodology, allows a much straightforward and complete study of systems with several non-connected foundations, and allows incorporating the flexibility of the foundation in the kinematic problem or analysing the soil response due to the vibration of the superstructure. The starting point of this proposal is the pile-soil interaction model presented in Álamo et al. (2016) where advanced Green's function were used to model layered soils and build an efficient model for the dynamic analysis of pile foundations. That paper also validated the model and the code for the problem of computing impedance functions. In the present paper, this idea has been completed and expanded to the seismic structure-foundation-soil problem by integrating the presence of a generic structure (modelled by beam and shell finite elements) and incident planar seismic waves with any angle of incidence. The proposed model is first satisfactorily validated against a previous BE-FE code for a base problem, highlighting the effects of considering the flexibility of the slabs on the structural response. Then, an application example of a pile-supported building structure founded on an inhomogeneous soil deposit on an elastic bedrock is considered. Besides the elements included in the application example, it is worth noting that the model can also represent, in its present state, additional features such as shear walls, inclined piles and any number and configuration of buildings, piles and foundations. The application example illustrates how the soil profile or the angle of incidence can affect significantly variables such as maximum efforts at pile heads, while other results, such as inter-storey drifts are dominated by the structural configuration.