Soil-Structure Interaction of Pile Groups for High-Speed Railway Bridges

STOCKHOLM SWEDEN 2018,

Soil-Structure Interaction of Pile Groups for High-Speed Railway Bridges

JOHANNES SEVERIN TOMMY STRAND

KTH ROYAL INSTITUTE OF TECHNOLOGY

In this thesis, the effects of soil-structure interaction on the dynamic properties of railway bridges are examined. Specifically, single-track bridges founded on pile groups and intended for high-speed traffic are studied. Due to the high speeds, these bridges must undergo a special dynamic analysis to ensure comfort and traffic safety. According to current practice, this analysis is carried out assuming fixed supports and any flexibility of the bridge piers is neglected.

The effects of flexible supports are in this work examined by representing the interaction between the pile foundation and the surrounding soil as a system of springs and dashpots. This system is described using impedance functions calculated from a finite element solution in frequency domain. The impedance functions are studied for load frequencies between 0 and 35 Hz and are calculated for six real pile groups with varying configurations and with three different soil types.

The impedance functions are then applied to simple two-dimensional bridge models in a complex eigenvalue study, where modal damping ratios and new natural frequencies for the first bending mode are calculated.

The greatest effects of the soil-structure interaction are seen for short single span bridges. The modal damping ratio of the first bending mode amounts to 6 % for a few cases and the largest decrease in natural frequency is 9 % or 1 Hz. However, for most cases, the effects are negligible.

Common to the largest observed effects is that they are given by pile groups of great depths. The soil modulus does not appear to be a decisive factor for either the damping or natural frequency, at least not for bridges with more than one span.

Keywords: Dynamics, SSI, Pile group, Foundation, HSR, Bridge, Impedance function, FEM, Damping, Complex eigenvalue analysis

I detta examensarbete undersöks hur jord-strukturinteraktion kan påverka de dynamiska egenska- perna hos järnvägsbroar. Speciellt studeras pålgrundlagda enkelspårsbroar avsedda för höghastig- hetstrafik. Till följd av de höga hastigheterna måste dessa broar genomgå en särskild dynamisk analys för att säkerställa komfort och trafiksäkerhet. Enligt nuvarande praxis sker denna ana- lys under antagande om fixa upplag och ingen hänsyn tas därmed till eventuell flexibilitet vid stöden.

Effekterna av eftergivliga upplag studeras här genom att jordens samverkan med pålgrundläggnin- gen betraktas som ett system av fjädrar och dämpare. Detta system beskrivs med hjälp av imped- ansfunktioner som beräknas ur en finita element-lösning i frekvensdomän. Impedansfunktionerna studeras för lastfrekvenser mellan 0 och 35 Hz och tas fram för sex verkliga pålgrupper med var- ierande konfiguration och med tre olika jordtyper. Impedansfunktionerna appliceras därefter på enkla tvådimensionella bromodeller i en komplex egenvärdesstudie, där modala dämpkvoter och nya resonansfrekvenser för första böjmoden beräknas.

De största effekterna av jord-strukturinteraktionen ses för korta enspannsbroar. Den modala dämpkvoten för första böjmoden uppgår till 6 % för några enstaka fall och resonansfrekvensen minskar med som mest 9 % eller 1 Hz. För de flesta fall är dock effekterna försumbara. Gemensamt för de största utslagen är att de ges av pålgrupper med stora djup. Jordens styvhet verkar inte vara en avgörande faktor för vare sig dämpningen eller resonansfrekvensen, åtminstone inte för broar med fler än ett spann.

Sökord: dynamik, jord-strukturinteraktion, pålgrupp, grundläggning, höghastighetståg, bro, im- pedansfunktion, FEM, dämpning, komplex egenvärdesstudie

This master’s thesis was initiated by the Department of Civil and Architectural Engineering at the Royal Institute of Technology, KTH, in conjunction with Tyréns AB.

First of all, we would like to thank our supervisor Ph.D. student Johan Lind Östlund for presenting this very interesting topic for us. He has shown a great commitment to our work and has been helpful in giving explanations and instructions throughout the process.

We would like to thank our examiner Prof. Costin Pacoste and Ph.D. student Abbas Zangeneh for the tutorial meetings we have had at KTH, where we have received meaningful input and guidance.

We would also like to thank the Bridge Department at Tyréns for giving us the opportunity to carry out this project at their office and for providing us with vital resources. Especially we want to thank the Head of Department Anna Jacobson and Ph.D. Mahir Ülker-Kaustell for their kindness and helpfulness.

Finally, we would like to thank our fellow thesis writers at the department, Arian Abedin, Wolmir Ligai, Christofer Schmied, Mattias Sedin and Benoît Van Gompel, for their support and great companionship during this work.

Stockholm, June 2018

Johannes Severin & Tommy Strand

Abstract iii

Sammanfattning v

Preface vii

Contents ix

1 Introduction 1

1.1 Background . . . 1

1.2 Aim and scope . . . 3

1.3 Limitations . . . 3

2 Theory 5 2.1 Dynamic properties of soil . . . 5

2.2 Model construction . . . 9

2.3 Complex eigenvalues and eigenvectors . . . 18

3 Method 21 3.1 Outline . . . 21

3.2 Description of the FE model . . . 22

3.3 Verification study . . . 24

3.4 Parametric study . . . 34

3.5 Complex eigenvalue analysis . . . 40

4 Results 45 4.1 Theoretical beam bridge . . . 46

4.2 Theoretical slab bridge . . . 48

4.3 Theoretical box bridge . . . 50

4.4 Real bridges . . . 50

5 Conclusions 51 5.1 Conclusions . . . 51

5.2 Discussion . . . 52

5.3 Future research . . . 54

Bibliography 55 Appendix 57 A Convergence analysis 57 A.1 Zheng et al. (2014) . . . 57

A.2 Padrón et al. (2012) . . . 59

B Verification study 60

B.1 Padrón et al. (2012) 2 × 2 pile group . . . 60

B.2 Padrón et al. (2012) 3 × 3 pile group . . . 62

C Receptance functions and impedance functions 64 C.1 Receptance functions for all DOF:s . . . 64

C.2 Vertical receptance functions for all pile groups . . . 67

C.3 Impedance functions for all DOF:s . . . 68

C.4 Vertical impedance functions for all pile groups . . . 70

D Results from the complex eigenvalue analysis 72 D.1 Theoretical beam bridges . . . 73

D.2 Theoretical box bridges . . . 77

D.3 Theoretical slab bridges . . . 81

D.4 Real bridges . . . 85

Introduction

1.1 Background

The Swedish Transport Administration (Trafikverket) is planning for a new high-speed railway network with the possibility of connecting the three largest cities in Sweden. The routes are shown in figure 1.1. The main purpose of the project is to increase the capacity of the railway system and to reduce the traveling time between the different regions. If a maximum allowed speed of 320 km/h is chosen, durations of 2 hours for the journey Stockholm-Göteborg and 2.5 hours for Stockholm-Malmö will be possible.

Figure 1.1: The planned Swedish high-speed railway network [1].

The construction will be made in stages, where Ostlänken between Järna and Linköping is the first. It comprises 150 km double-track railway and around 200 bridges. The aim is for Ostlänken to be completed in 2028 and the entire high-speed railway network in 2035 [1].

1.1.1 Dynamic analysis on railway bridges

Bridges on high-speed railways are exposed to high dynamic loads from the running trains. If the maximum allowed speed is to exceed 200 km/h, an analysis of the dynamic effects is required [2]. The relevant verifications are stated in the Eurocodes (SS-EN 1991-2 6.4) and include, among others, bridge deck accelerations and displacements [3].

The Eurocodes suggest values for the structural damping that may be used in the dynamic ana- lysis. The appropriate value is to be selected depending on the type and span length of the bridge in question, according to the table 1.1. No distinction is made between different sources of damping.

Table 1.1: Damping ratios, ξ (%), to be used in dynamic analyzes [3].

Bridge type Span L < 20 m Span L ≥ 20 m

Steel and composite ξ= 0.5 − 0.125(20 − L) ξ= 0.5 Prestressed concrete ξ= 1.0 + 0.07(20 − L) ξ= 1.0 Filler beam and reinforced concrete ξ= 1.5 + 0.07(20 − L) ξ= 1.5

1.1.2 Soil-structure interaction and support flexibility

Common practice in bridge design is to assume the deck to be supported on fixed supports. For structures founded on soil, there is however a certain support flexibility which will affect the dynamic properties of the bridge. The effects of soil-structure interaction (SSI) are difficult to estimate and analyzes containing the entire system of bridge, foundation and soil are computa- tionally costly.

An approach to reducing the computational time is to disconnect the soil and foundation from the bridge model and to calculate the impedance functions of the supports in a separate model. The impedance functions represent the SSI and consist of expressions for the dynamic stiffness and damping of the soil-foundation system. This approach reduces the size of the model and allows for the impedance functions to be used on other bridge models as well. A comparison between this approach and the assumption of fixed supports is shown schematically in figures 1.2 and 1.3.

Springs and dashpots are used to represent stiffness and damping, respectively.

Figure 1.2: Bridge model with fixed supports.

+ =

Figure 1.3: The impedance functions representing SSI are being attached to a bridge model. The result is a bridge model with flexible supports.

1.1.3 Pile group foundations

In a previous study from KTH [4] the effects of support flexibility on the dynamic performance of railway bridges were examined for two types of foundations. The results indicated that the effects of SSI might be significant for bridges founded on pile groups. It is therefore of interest to further examine the factors contributing to the interaction between the soil and piles.

1.2 Aim and scope

The aim of this thesis is to investigate the SSI of pile groups and to study the effects on the dynamic performance of railway bridges. The process can be divided into steps:

• A verification study, in which the impedance functions are compared to analytical and nu- merical solutions from the literature.

• A parametric study, in which the influence on the impedance functions from different soil parameters and pile group configurations is examined.

• A complex eigenvalue analysis, in which the effects of the impedance functions will be de- termined in terms of modal damping ratio and change of natural frequency.

1.3 Limitations

A list of the most important limitations and assumptions are presented below.

• The soil is treated as a homogeneous, isotropic, viscoelastic single layer medium overlying an infinitely stiff bedrock.

• The piles are end-bearing and assumed fully fixed to both the underlying rock and to the pile cap.

• Fully bonded contact is assumed between the piles and the surrounding soil.

• No contact is assumed between soil and pile cap.

• A frequency-independent structural damping model is applied for the soil medium.

• The dynamic properties of the bridges are calculated only for the first vertical bending mode.

Theory

2.1 Dynamic properties of soil

In this thesis, the soil is considered as semi-infinite since it is bounded in the vertical direction by the surface and the bedrock, and unbounded in the horizontal direction.

2.1.1 Basic soil concepts

The equations regarding the basic soil concepts are thoroughly explained in [5]. The shear modulus Gfor a material is the ratio of shear stress τ and angular displacement γ when subjected to shear.

For a linear elastic material the equation is:

G= τ

γ [Pa] (2.1)

However, for a soil, the shear modulus decreases with increasing deformation and the relationship between stress and deformation is hyperbolic. The elastic modulus E is the ratio of axial stress and deformation:

E=σ_x

ε_x [Pa] (2.2)

Poisson’s ratio is the negative ratio between transverse strain and axial strain:

ν= −ε_z

ε_x [−] (2.3)

If the elastic modulus and Poisson’s ratio of a soil is known, the shear modulus of a soil is:

G= E

2(1 + ν) [Pa] (2.4)

When calculating soil settlements due to load on the ground surface, the compression modulus M is commonly used. It is determined by increasing the effective vertical pressure in a sample while preventing horizontal deformations:

M = E(1 − ν)

(1 + ν)(1 − 2ν) [Pa] (2.5)

2.1.2 Waves in an elastic semi-infinite soil

The theory in section 2.1.2 is based on writings by Bodare [6], unless otherwise cited. Inside a soil volume two types of waves can exist, P-waves and S-waves , with different propagation velocities.

These waves exist separately and are independent of each other.

The P-wave propagates by dilation (change of volume) where the material particles move parallel to the direction of the wave propagation. This occurs through compression and rarefacton in the material. P stands for ’primary’ since the P-wave travels fastest in the ground and is first wave to be recorded at an observation point. The propagation velocity of a P-wave is:

c_p= s

M

ρ. [m/s] (2.6)

S-waves propagate through shear deformation, where the material particles move perpendicular to the propagation direction. S stands for ’secondary’ because it propagates half as fast as the P-waves. The S-wave cannot propagate in water, which implies that the velocity is the same under the ground water surface. Since it propagates by shear deformation it is dependent on the shear modulus G. The velocity of the S-wave is hence determined by:

c_s= s

G

ρ [m/s] (2.7)

λ λ

Figure 2.1: The P-wave transmits energy by change of volume, and the S-wave by change of shape.

According to [7] typical S-wave velocities are 0-250 m/s for clays, 200-400 m/s for sand and gravel and 200-700 m/s for moraine.

When P-waves and S-waves are reflected at the ground surface so-called Rayleigh waves, or R- waves, are generated. R-waves resemble the surface waves that are created when a rock is thrown into water [7]. The wave propagation velocity is approximately 90 % of the S-wave’s and they are hence difficult to separate. The velocity is calculated:

c_R= 0.87 + 1.2ν

1 + ν · c_S [m/s] (2.8)

The soil is typically layered with varying properties, due to sedimentary depositions. Even a seemingly completely homogeneous soil in fact becomes inhomogeneous, as the pressure from the overlying soil causes the shear modulus to increase with the depth. Reflection and refraction occurs when a wave enters a new layer. Love waves, or L-waves occur when S-waves hit the interface of layer with a significantly higher shear modulus G (e.g. rock). If the angle of incidence is large enough, a total reflection of the waves will be obtained and all energy remains in the interface.

L-waves can therefore propagate long distances.

Figure 2.2: Typical R-wave and L-wave.

Frequency domain

The characteristics of a wave can be described in both time domain and in frequency domain.

In the time domain it is possible to evaluate the maximum and minimum values for the soil displacements, while frequency domain is used to evaluate the dominant frequencies [7]. Basic knowledge of transition from time domain to frequency domain using Fourier transformation is assumed, see figure 2.3.

A₁

A₂

t U(t)

f A(f)

f₀ 6f₀

A2

A₁ Fourier

Transformation

Figure 2.3: Fourier transformation from time domain to frequency domain of a complex signal comprising two different frequencies with different rates and amplitudes.

2.1.3 Damping properties of the soil

All waves in the soil are damped, due to wave propagation and friction. The damping depends on the soil properties and frequency. In a soil there are two types of damping that contribute to the decay of a dynamic process: material damping and radiation damping. The theory in section 2.1.3 is based on writings by Bodare [6], unless otherwise cited.

Material damping

For small shear strains in a soil (< 0.001 %) the material is linearly elastic and its ability to dissipate energy is negligible. Hence, there is no active material damping in the soil. However, at moderate strains (0.01 − 1 %) most materials exhibit both elastic and plastic properties (elastoplasticity).

A residual deformation is obtained after unloading. Increased strain implies decreasing shear modulus G and increasing damping. The damping is generated from friction between the soil particles and is independent of the frequency. When plotting the cyclic shear strain γ against the shear stress τ a loop is obtained that encloses a surface, see figure 2.4. The surface is called the hysteresis surfaceand its area corresponds to dissipated energy ∆W for each cycle. The maximum potential energy W corresponds to the triangular area under the straight line. The damping ratio ξis the ratio of dissipated energy and potential energy:

ξ= 1 4π ·∆W

W [−] (2.9)

γ

τ Gmax G

τa

γa

ΔW W Initial

deformation curve

Figure 2.4: Hysteresis loop showing the areas of ∆W and W . τa and γa is the maximum amplitude of shear tress and shear strain, respectively. The small degradation of the shear modulus is shown by the shift from initial shear modulus Gmaxto G.

For even higher strains (> 1) % a degradation of the cyclic stiffness arises. That implies that the shear modulus G is reduced for each cycle which makes the soil much more complicated to model and calculate.

Radiation damping

Radiation damping is caused by the waves spreading over a larger volume, since the total energy over the expanded wavefront must be constant [7]. The wave amplitude decreases with increasing distances from the source according to:

A₂ A₁ =r₂

r₁

⁻ⁿ

e^{−α(r2−r1)} [−] (2.10)

where n = 0.5 for L-waves and R-waves, n = 1 for P-waves and S-waves. r1 is the inner radius and r2 the outer, respectively. α is the material damping factor and is dependent on the damping ratio ξ, frequency f, and S-wave propagation velocity cS:

α= 2πξf

c_S [m⁻¹] (2.11)

This shows that the radiation damping is frequency dependent, unlike material damping. Higher frequencies are more damped than lower frequencies.

2.2 Model construction

The response of soils exposed to dynamic loads is not easily computed. The cyclic behavior under random or cyclic loading is important to model. In this section, some measures of simplifications and assumptions are introduced in order to facilitate the analysis.

2.2.1 Linear viscoelastic model

As long as the cyclic shear strains are small, 0.01 − 1 %, the properties of the soil can be described by a constitutive model using classical theory of linear viscoelasticity, called linear viscoelastic model. Viscoelasticity means that a material exhibits both viscous and elastic properties during deformation. Viscosity is the linear resistance between shear stress G and strain γ. Linear elasticity means that the material returns to its original position after the load has been removed. Hence, in a viscoelastic model, the shear-strain relationship is linear, but also damped. It is proven that there is always a certain degree of damping in soils. Thus, this model can be used to represent the properties of the soil, even when the soil exhibits non-linear behavior, where otherwise the damping plays a major role. The theory and derivation of loss coefficient in section 2.2.1 is collected from Ishihara [8].

Loss coefficient

A commonly used soil parameter is the loss coefficient η, which indicates how much of the vibration energy that is converted to heat through friction. It is dependent on both the elasticity and the energy dissipating properties of the soil. The derivation of the loss coefficient is as follows. A harmonic shear stress applied to a viscoelastic body is denoted:

τ = τasin (ωt) [Pa] (2.12)

where τ_a is the amplitude, t is the time, and ω is the angular frequency. The shear strain will respond with the same frequency, but with a small time delay δ:

γ= γasin (ωt − δ) [−] (2.13)

where γ_texta is the strain amplitude. The strain response, which is commonly described by the ratio τa/γ_a, now has to consider the time delay δ. A method to solve this is to introduce complex variables. The stress and strain can then be presented:

τ_R= τ_acos (ωt) and γ_R= γa cos (ωt − δ) (2.14) respectively. Index R indicates the real value in a complex number. Thus, the complex represent- ation of the stress and strain can be expressed:

τ = τR+ iτ = τasin (ωt) + iτacos (ωt) [Pa] (2.15) and

γ= γR+ iγ = γasin (ωt − δ) + iγacos (ωt − δ) [−] (2.16) respectively. By using Euler’s formula a more slender formulation is obtained:

τ= τ_ae^iωt and γ = γ_ae^i(ωt−δ) (2.17)

The strain response τ_a/γ_a is now described as:

τ γ = τ_a

γ_ae^iδ= τ_a

γ_a(cos δ + i sin δ) [Pa] (2.18) The parameters elastic modulus µ and loss modulus µ⁰ are now introduced. The elastic modulus represents the elastic or instantaneous response and the loss modulus represents the dissipated energy of a viscoelastic body. They are obtained by admitting:

µ= τ_a

γ_acos δ and µ⁰= τ_a

γ_asin δ (2.19)

μ'γa μγa

γ τ

ΔW W

Figure 2.5: Hysteresis loop for a viscoelastic model. A simple way to graphically determine the loss coefficient η is to read off the values of the elastic modulus µ and loss modulus µ⁰.

The total strain response τ_a/γ_aconsidering both the elastic and damping properties can be presen- ted by the complex modulus, denoted:

µ^∗= µ + iµ⁰= τ_a

γ_a [−] (2.20)

alternatively expressed:

|µ^∗|=p

µ²+ µ⁰²= τ_a

γ_a [−] (2.21)

The loss coefficient η is the ratio of the elastic modulus and loss modulus:

η= µ⁰

µ = tan δ [−] (2.22)

and δ is the time delay.

Frequency-independent Kelvin-Voigt model

The viscoelastic behavior of a material is visualized by simple models containing springs and dashpots. The springs represent the elastic properties and the dashpots represent the damping properties. They can be connected to each other either in parallel or in series. A viscoelastic model in which the spring and dashpots are connected in parallel is in the literature referred to a Kelvin-Voigt model. The damping is usually treated as frequency-dependent. By introducing the spring constant Gand dashpot constant G⁰ the total stress for a Kelvin-Voigt model is:

τ= (G + iωG⁰)γ [Pa] (2.23)

Compared to eq. 2.19 the elastic modulus µ and the loss modulus µ⁰ are simply:

µ= G and µ⁰= G⁰ω (2.24)

The loss factor η compared to eq. 2.22 is now:

η= G⁰ω

G = tan δ [−] (2.25)

However, we know that the material damping in soils are nevertheless considered independent of the frequency. Therefore, the application of a pure Kelvin-Voigt model may in practice only be useful for special cases with low frequency or negligible damping. In order to take advantage of the benefits of the Kelvin-Voigt model, it is necessary to introduce a special type of damper that is independent of the frequency. The easiest way to meet the requirements is to use the frequency-independent so-called non-viscous type Kelvin-Voigt model, see figure 2.6. The total stress is obtained simply by removing the frequency factor ω in equation 2.23, so that:

τ= (G + iG⁰)γ [Pa] (2.26)

The model now consists of a spring and a frequency-independent dashpot that also are connected in parallel. The equation 2.26 can be interpreted as a stress-strain ratio where the stress consists of two parts: one where the stress coincides with the strain, τ₁ = Gγ, and the other with a phase shift of 90^◦, τ₂ = G⁰γ. The stress-strain ratio contains an imaginary term, which has no physical implication. Compared to equation 2.24 the elastic modulus µ and the loss modulus µ⁰ corresponding to a non-viscous type Kelvin model are now:

µ= G and µ⁰= G⁰ (2.27)

which, compared to equation 2.25, give frequency-independent loss factor:

η= G⁰

G = tan δ [−] (2.28)

τ

τ₁ τ₂

γ

G G'

F(ω)

Figure 2.6: A model showing a viscoelastic medium containing a spring and a frequency-independent dashpot, illustrating the properties of a frequency-independent Kelvin-Voigt model. τ1is the stress carried by the spring and τ2is carried by the dashpot.

Structural damping

It is considered more advantageous to create models with a linear damping behavior, rather than approximate the viscoelastic damping properties of an actual soil. According to Chopra [9] it is possible to construct a damping coefficient that is equal to the combined effects of all damping contributions. The equivalent damping coefficient for viscous damping ξ_eq(also denoted structural damping) is based on actual measured responses of complex systems with many degrees of freedom (DOF:s) exposed to a harmonic load oscillating in the same frequency as the natural frequency of the structure. For those special cases the structural damping is related to the loss coefficient in the simple manner:

ξ_eq= η

2 [−] (2.29)

The structural damping ratio is not accurate for other frequencies besides the natural frequency of the structure but is still a satisfactory approximation. As seen in figure 2.7 the hysteresis loop takes an arbitrary shape, due to real-life force-displacement relationship.

γ τ

ΔWeq

W_eq

Figure 2.7: Hysteresis loop of actual force-displacement relationship of a harmonic excitation, determined from experiment. ∆Weqis the measured cyclic dissipated energy of an actual structure, and Weq= ku²₀/2, k being the stiffness determined by experimentation and u0 being the displacement amplitude.

2.2.2 Soil-structure interaction

Soil-structure interaction (SSI) arises in a system consisting of two parts with significantly different damping properties. The stiffer structure, or the more flexible soil, the more the SSI increases.

For static loads, SSI is not taken into account. But as the load frequency increases, dynamic conditions arise, and SSI should normally be taken into account. The theory in section 2.2.2 is based on writings by Wolf [10].

The dimensions of a structure are finite, which implies that it is relatively straightforward to generate a dynamic finite element model that comprises a finite number of nodes with different DOF:s. For higher frequencies, SSI is taken into consideration which means that also a part of the soil region must be modeled. The soil, however, is a semi-infinite medium. For a static load, one may introduce a fictitious boundary at a sufficient distance where the response is expected to die out. In this way, one obtains a finite domain, comprising the soil and the structure, that can be calculated in a similar way as the structure. However, for a dynamic load, the fictitious boundary will reflect waves, rather than allowing them to pass through and propagate to infinity.

Furthermore, creating a model where the structure is embedded or integrated in the ground leads to many DOF:s, which requires a lot of computational work.

It is, consequently, preferable to analyze the soil as a dynamic subsystem. The two subsystems, the soil and the structure, can then be analyzed separately using two different methods. The dynamic stiffness is frequency-dependent and complex. This implies that linear SSI analysis is best treated in the frequency domain. The benefits of operating in the frequency domain are important. By means of so-called impedance functions, the natural frequencies, with the associated damping ratio, can be determined and controlled. The method of separating the soil from the structure breaks down the complicated soil-structure system into more manageable parts that are easier to control, and the impedance function of the soil does not need to be recalculated for other load cases. The same applies if the structure is modified.

2.2.3 Impedance function

The derivation of and mathematical relationship of an impedance function in 2.2.3 is from Gazetas [11]. Impedance, written in the time domain, is the ratio of force F (t) and displacement u(t), at the nodes that are in direct contact with the structure:

Z(t) = F(t)

u(t) [N/m] (2.30)

As the impedance is frequency-dependent it is better described in the frequency domain:

Z(ω) = F(ω)

u(ω) [N/m] (2.31)

The impedance function can also be represented by a complex formulation:

Z(ω) = k(ω) + iωc(ω) [N/m] (2.32)

The function k(ω) (N/m) is the dynamic stiffness and is also described as the real part of the impedance function Re(Z). It reflects a combination of the static stiffness and the inertia of the soil. For low frequencies the dynamic stiffness decreases parabolically, but increases again at a certain frequency, due to resonance.

The function ωc(ω) (N/m) is the product of the frequency ω and the dashpot coefficient c(ω) (Ns/m) and is described as the imaginary part of the impedance function Im(Z). For low fre- quencies the imaginary part is low, but increases for higher frequencies, due to resonance. The damping coefficient reflects the total damping generated by both radiation damping and material damping.

F(ω)

k(ω) c(ω)

u(ω)= k(ω)+iωc(ω)F(ω)

Figure 2.8: Illustration of the impedance of the soil, using a Kelvin-Voigt model. The displacement, u(ω), is the ratio of the force and impedance.

For a unit load F (ω) = 1 N the impedance function is obtained by the inverse of displacement u(ω) = 1/Z(ω) (m), commonly referred to as the receptance. In the frequency domain the recept- ance is described in so-called frequency response functions (FRF:s).

Derivation of impedance function from EOM

The equation 2.32 is applicable for all foundation-ground systems. The interpretation of k(ω) and c(ω) as dynamic stiffness and damping coefficients, respectively, can be justified by derivation of impedance using the equation of motion (EOM), which in the time domain is known as:

M¨u(t) + C ˙u(t) + Ku(t) = F (t) [N] (2.33) where ¨u(t) is the acceleration ˙u(t) is the velocity and u(t) is the displacement. M, C and K are here treated as constant matrices. By introducing the Euler’s identity and describing the displacement u in frequency-dependent amplitudes for a constant frequency ω one obtains:

u(t) = u(ω)e^iωt [m] (2.34)

The derivative and the second derivative is:

˙u(t) = iω · u(ω)e^iωt and ¨u(t) = −ω²· u(ω)e^iωt (2.35) respectively. Now the EOM in the frequency domain may be written as:

M ·(−ω²) · u(ω)e^iωt+ C · iω · u(ω)e^iωt+ K · u(ω)e^iωt= F (t)e^iωt [N] (2.36) After simplification one obtains the equation:

(−ω²M + iωC + K)u(ω) = F (ω) [N] (2.37) Hence, the impedance function is:

Z(ω) = (−ω²M+ iωC + K) = F(ω)

u(ω) [N/m] (2.38)

2.2.4 Finite element method for solids

The purpose of section 2.2.4 is to give an overall explanation of the most important notions of FEM for solid materials with isoparametric elements. Isoparametric means that an element face is not necessarily rectangular. Most of the material and equations presented in section 2.2.4 is based on Chopra [9] and Felippa [12] The more specific cases regarding infinite elements and embedded beams, are from [13].

With FEM it is possible to formulate the EOM for complicated structures. An arbitrary solid structure is discretized in many small finite elements that are connected to each other in nodes. A node can have several DOF:s, both translational and rotational, depending on geometry. For an isoparametric solid element, however, there are only three DOF:s, namely translational displace- ment in the x, y and z direction.

Derivation of the EOM

For a two-dimensional beam element, a shape function describes how the beam will be deformed between two nodes. For those cases, it is easy to illustrate shape functions, but for three- dimensional cases the shape functions are based on assumed relationships between displacements in the element’s internal points and displacements in the nodes. This makes it easy to handle, but introduces some approximations. The errors can be reduced by reducing the element size or increasing the number of finite elements. A vector consisting of shape functions for each node, i, of a three-dimensional element with n number of nodes is created:

N=





 N₁

...

N_n





 (2.39)

For each element, the element stiffness matrix K_l, the element mass matrix M₁, and the ele- ment force vector F_l are generated. Index l indicates that the matrices are described with local coordinates. The vector u represents all node displacements for each DOF:s in one element:

u= ux1 u_y1 u_z1 . . . u_xn u_yn u_zn^T (2.40) A so-called B matrix is created comprising different derivatives of all shape functions Ni with respect to x, y and z:

B =



























δN₁

δx 0 0 . . . ^δN_δxⁿ 0 0 0 ^δN_δy¹ 0 . . . 0 ^δN_δyⁿ 0 0 0 ^δN_δz¹ . . . 0 0 ^δN_δzⁿ

δN₁ δy

δN₁

δx 0 . . . ^δN_δy^{n δN}_δxⁿ 0 0 ^δN_δz^{1 δN}_δy¹ . . . 0 ^δN_δz^{n δN}_δyⁿ

δN1

δz 0 ^δN_δx¹ . . . ^δN_δzⁿ 0 ^δN_δxⁿ



























(2.41)

The strain for each node is calculated ε = Bu.

For an isotropic material, the relationship between the stress and strain is constant and therefore the matrix of the elastic modulus can be described:

E=























Eˆ(1 − ν) Eνˆ Eνˆ 0 0 0 Eνˆ Eˆ(1 − ν) Eνˆ 0 0 0 Eνˆ Eνˆ Eˆ(1 − ν) 0 0 0

0 0 0 G 0 0

0 0 0 0 G 0

0 0 0 0 0 G























(2.42)

where ˆE= (1−2ν)(1+ν)^E . The local element stiffness matrix is then:

K_l =Z

Vl

B^TEB dV (2.43)

where V_l is the volume of the element. To derive the so-called consistent mass matrix, the shape functions are rearranged to:

N=





N₁ 0 0 . . . Nn 0 0 0 N₁ 0 . . . 0 Nn 0 0 0 N₁ . . . 0 0 Nn



 (2.44)

so that the consistent element mass matrix is:

M_l=Z

V_l

ρN N^T dV (2.45)

where ρ is the density of the material. The consistent mass matrix is not diagonal, which com- plicates the computation process. It is simplified by assuming that the distributed mass in the element can be lumped together as point masses at the nodes. This matrix is called lumped-mass matrix. The damping matrix can not be calculated based on the dimensions of the structure, member sizes or the damping of the structural material. Even though the damping properties of each member are known, energy dissipation is not taken into account. The damping matrix can be determined based on available data for similar structures, or table cases. If the external forces are applied along the three DOF:s at the nodes of a finite element, the force vector of an element can be written directly. All local matrices that apply to an element are assembled into proper locations in global matrices that apply to the entire system K_g, M_g, C_g and F_g, respectively.

Finally, the equation of motion is formulated:

M_g¨u + C_g˙u + K_gu= F_g(t) (2.46)

Some element properties

Three-dimensional finite elements have three standard geometries: tetrahedron, pentahedron and hexahedron, see figure 2.9. In addition there are two nonstandard geometries used to fill regions that are built up of hexahedron elements. Elements with both standard and nonstandard geomet- ries can be refined with additional nodes, preferably between each corner, to get a more accurate stress analysis. Unlike linear formulation, one obtains a quadratic formulation of the displace- ments in each node, see figure 2.10. Quadratic formulations provide more accurate results for larger element sizes but requires more computational time.

Figure 2.9: Standard shapes for three-dimensional elements, with linear formulation. In order: tetrahed- ron, pentahedron and hexahedron.

Figure 2.10: Standard shaped three-dimensional elements with quadratic formulation.

In order to obtain good representation of the wave movements in a FE model, the smallest occur- ring wavelength should be possible to be described by at least 10-12 elements (i.e. 10-12 nodes) for linear formulation, or 4-5 elements (or 8-10 nodes) for quadratic formulation. The maximum size, depending on element geometry, is for linear elements approximately ^λ₈ and for quadratic elements ^λ₄, respectively [14].

Infinite elements are useful in infinite or semi-infinite domains, for instance in cases with wave propagation in a soil, where the waves propagate towards infinity. An outer boundary of a model, consisting of regular finite elements, will reflect waves back into the model, resulting in super- position of the transmitted and reflected waves. One solution would be to construct a boundary, consisting of finite elements, far away from the region of interest. This, however, entails a lot of elements and a long computational time. Therefore, the use of infinite elements is an effective method of dynamic analysis. Boundaries with infinite elements give best results for cases where the dominant direction of the S-waves are orthogonal to the boundary. Therefore, the infinite boundary should be placed at a suitable distance from the region of interest so that the angle of incidence becomes as orthogonal as possible, and the amplitude small enough.

S-wave S-wave S-wave

Reflection Reflection

Finite boundary

Finite boundary Infinite boundary

<90° ≈90°

1. 2. 3.

Figure 2.11: Three approaches to model the boundary of a semi-infinite finite element model with wave propagation.

The first case in figure 2.11 is the worst, where the radius is small which leads to deficient reflection of the waves, due to small angle of incidence. The second case has a satisfactory angle of incidence but still produces reflection and requires many elements, hence longer computational time. The third and best case has an infinite element boundary, which requires less amount of elements and produces no reflection.

In some cases, it is desirable to model a two-dimensional beam of any cross-sectional dimension, length and property, inside a three-dimensional solid material. An embedded beam in a solid is a two-dimensional beam divided into elements with connecting nodes embedded in a group of so- called host elements belonging to a solid body. Each node in the beam take over the translational degrees of freedom of the solid host element surrounding that node. Embedded elements can keep their own rotational DOF:s, and these rotations are not locked to the host element. One beam can have several nodes within a host element.

Embedded beam

Host element

Host element nodes

Embedded nodes coupled to host

Figure 2.12: Illustration of a two-dimensional embedded beam in solid host elements.

2.3 Complex eigenvalues and eigenvectors

Eigenvectors determine the natural mode shapes of vibration by defining the displacement of each node in a multi-degree of freedom system vibrating at its natural frequency. For undamped systems the eigenvector is usually denoted Φn, and for damped systems it is denoted by complex-conjugate pairs, Ψn and Ψn. Index n stands for the specific mode. The complex conjugate of a complex number is a number with equal real part and imaginary part and has the same magnitude, but opposite sign and is marked with an overline. The theory in section 2.3 is from Chopra [9], unless otherwise cited. For a system with N degrees of freedoms the equation of motion for free vibration and damping is:

M¨u + C ˙u + Ku = 0 (2.47)

To solve the equation following is admitted:

u(t) = Ψe^λt (2.48)

where Ψ is the eigenvector and λ is the eigenvalue, associated with the natural frequency. Equation 2.48 in equation 2.47 yields the so-called complex eigenvalue problem, which has no real roots:

(λ²M+ λC + K)Ψ = 0 (2.49)

The eigenvector Ψ is best solved by reducing the system of N second-order differential equations into a system of 2N first-order differential equations. For cases with low damping the roots will occur in complex-conjugate pairs, where the real part is either negative or zero. For systems with N degrees of freedom there will be N pairs of eigenvalues. The solution to the equation 2.49 gives the eigenvalues:

λ_n= −ξnω_n+ iωnD and λn= −ξnω_n− iωnD (2.50) where the natural frequency for damped systems is:

ωnD = ωn

p1 − ξn² (2.51)

The natural frequencies ωn and the modal damping ξn is related to the eigenvalues as fol- lows:

ωn = |λn| and ξn= −Re(λn)

|λn| (2.52)

For each pair of eigenvalues there is a corresponding complex-conjugate pair of eigenvectors, with a real part and an imaginary part:

Ψn = Φn+ iχn and Ψn= Φn− iχn (2.53) where Φn and χn are real valued vectors with N elements each. Φn is the eigenvector for an undamped system and χn is a vector representing the damping contribution. Hence, for χn = 0, Ψn = Ψn = Φn.

Modal assurance criterion(MAC) is a useful tool in complex eigenvalue analysis [15]. It is a method based on the method of least squares and is used to investigate statistical relationships between two related mode shapes for a structure exposed to many frequencies. One example of related mode shapes is where one mode shape is provided analytically and the other one experimentally. Another example is a mode shape before and after a change in the structure caused by modification, e.g.

added impedance functions. MAC is most sensitive to large differences and relatively insensitive to small differences of mode shapes, which indicates statistical consistency between them. It is calculated as a normalized scalar product of two sets of compared eigenvectors, ΦAand ΦB:

M AC(r, q) = |{Φ_A}^T_r{Φ_B}_q|²

({ΦA}^T_r{ΦA}r)({ΦB}^T_q{ΦB}q) (2.54) in which q and r are modes. The form of coherence function is recognized, indicating the causal relationship between the eigenvectors. The MAC number is bounded between 0 and 1, where 1 indicates fully consistent mode shapes and a number close to 0 indicates that the mode shapes are not consistent.

Method

3.1 Outline

The soil-structure interaction was investigated by creating an FE model consisting of a pile group and the surrounding soil and by studying the steady-state response to a harmonic unit load applied to the pile group. From the obtained displacements, the corresponding impedance functions were calculated. The SSI effects on bridges were determined by applying the impedance functions to a set of beam models and by performing a complex eigenvalue analysis.

These steps are described in detail in the following sections, where section 3.2 deals with the features of the FE model, section 3.4 the computation of impedance functions and section 3.5 the complex eigenvalue analysis. A flowchart of the procedure is shown in Figure 3.1.

Abaqus

• Steady-state dynamics direct procedure Input

• Model properties

• Pile configuration

Frequency Response Functions Python

• Inverting

• Sorting

• Plotting Impedance

Functions

Matlab

• Simulation of bridges

• Comparing mode shapes

• Plotting

Results

• Damping ratios

• Natural frequency

Figure 3.1: Flowchart of the procedure for determining the SSI effects on railway bridges.

In terms of software, Abaqus was the main tool for modeling, simulation and computation of the response of the soil and piles to the applied dynamic loads. Abaqus (or Abaqus FEA) is a product suite by the software company Dassault Systèmes (3DS) for finite element analysis and computer- aided engineering. All products use the open-source scripting language Python for scripting and customization. For this work, Brigade Plus by Scanscot Technology was used, which is a software series based on Abaqus focused on structural, bridge and civil engineering. The computer programs and languages Matlab and Python was used for data processing and plotting.

3.2 Description of the FE model

The response of the soil volume and the pile group due to an applied harmonic unit load is studied in a steady-state dynamics direct procedure. A disc-shaped model is chosen with varying heights depending on the pile lengths. The radius must be large enough to prevent the so-called ripple effect, which gives an unstable description of the FRF:s and impedance functions. Therefore, different radii have been examined, resulting in a final radius which is 200 m. The pile group is placed in the center of the model, without the pile cap, and each pile extends to the bottom of the model. A reference point for load application and for information about the output is placed in the center of the pile group at the level of the ground surface. The FE model is divided into the four regions:

• Piles

• Soil (mid region)

• Soil (main region)

• Boundary

as indicated in figure 3.2. The mid region is the region that encloses all the piles and has approx- imately equal element sizes. Its radius also includes an additional meter from that pile toe which is located furthest away from the origin. In order to limit the number of elements and computational time, there is a radial dilution (so-called bias) of the element mesh in the main region starting from the mid region interface. The boundary region consists of infinite elements.

Boundary Piles

Mid region

Main region

Figure 3.2: Screenshot of the entire FE model indicating the four regions.

3.2.1 Element types

The entire soil region is modeled with tetrahedrons as they provide a good description of the wave motions in all directions. A convergence analysis showed that the use of quadratic formulation is preferable, in Abaqus denoted C3D10. The piles are modeled with the same element type as the adjacent soil region, but with other material properties. The boundary region is modeled with infinite elements with hexahedrons and quadratic formulation, CIN3D8.

3.2.2 Boundary conditions and constraints

The entire bottom surface of the soil region is prevented from moving in all directions.

Each pile head surface is coupled to each other with a so-called rigid body constraint. This implies that they move together as a rigid body, which they actually do when they are connected to a pile cap. This body containing the pile heads is free to move in all translational and rotational directions and is coupled to the reference point, see figure 3.3.

The piles extend all the way to the bottom of the model, and the pile toes are fixed.

Since the soil region and the boundary region are modeled with different element geometries, they are incompatible with each other. However, so-called tie constraints are created between the nodes of the soil elements and the infinite elements.

Reference point

Figure 3.3: Screenshot of the a pile group indicating the pile head surfaces which are tied to the reference point with a rigid body constraint.

3.2.3 Distribution of the element sizes

Assigned mesh size in the mid region is 0.5 m, and increases to 20 m at the edge of the model.

The piles are modeled with solid elements, with assigned mesh size 0.5 m (the same as the mid region). Due to the geometries of the piles the base pile elements have the side length of approx.

0.22 m. See section 3.3.5 for further details.

The maximum permissible element size in the mid region is λmin/5, as discussed at the end of section 2.2.4. The element size is governed by the pile dimensions and never exceed λ_min/5.

Reference point

Mid region Main region

Piles

Main region

Figure 3.4: Section of the mid region showing the piles in relation to the adjacent soil. Pile elements govern the size of the surrounding elements in the mid region. The bias in the main region is apparent.

3.2.4 Implications of a linear viscoelastic model

The soil volume is evaluated according to a frequency-independent linear viscoelastic Kelvin-Voigt model, which means that the soil is linearly elastic and damped and that the loss factor η is the ratio between the dashpot constant µ⁰ and the spring constant µ. The model is frequency-independent in the sense that the material damping is constant for all frequencies. Therefore, the structural damping ξ_eqis used in Abaqus, which is half the loss factor η. One alternative would be to use viscoelastic damping, but this requires considerably more knowledge about the soil parameters, which furthermore are frequency dependent. The damping ratio used in the parametric study is ξ_eq= 2 %.

3.3 Verification study

The Abaqus FE model was verified by comparing the obtained impedance functions to solutions from previous studies. The purpose was to ensure the validity of the results gained in the para- metric study. This verification study includes comparisons to the following:

• Kobori et al. (1971), presented in section 3.3.1

• Boundary element method, presented in section 3.3.2

• Zheng et al. (2014), presented in section 3.3.3

• Padrón et al. (2012), presented in section 3.3.4

Furthermore, the stiffness of the three-dimensional solid piles in the FE model was compared to standard table cases in section 3.3.5.

In general, the analyzed parameters in order to reach convergence were total radius, mid-region radius, element types, element sizes and bias. Comparisons between piles modeled with embedded beam elements and solid elements were also made.

3.3.1 Verification against Kobori et al. (1971)

Kobori et al. (1971) [16] is a semi-analytical solution of a 4 meter deep viscoelastic, homogeneous and isotropic soil on a rigid bottom, subjected to a vertical harmonic load evenly distributed over a flexible quadratic area of 4 m². The impedance functions are presented as functions of the non- dimensional frequency a₀= ωB/c_s, where B = 1 m. The real and imaginary part of the impedance function is presented with equivalent stiffness and damping coefficients such that:

K_ev= Re(u)

Re(u)²+ Im(u)² and C_ev= − Im(u)

a₀(Re(u)²+ Im(u)²) (3.1) and is independent of the height H and elastic modulus E. Used parameters for the verification model is:

H = 4 m E = 400 MPa ρ = 2000 kg/m³ ν = 0.1

η = 0.1

H E, ρ, ν, η H E, ρ, ν, η

F/m² F

Figure 3.5: Assumed conditions for verification according to Kobori et al. (1971) with a flexible plate.

This gives a shear wave velocity of c_s≈283 m/s. The structural damping ξ_eq= 5 % is used. The result was obtained for the radius of the mid region set to 4.5 m and mesh size 0.25 m, while the total radius was 200 m with an assigned mesh size of 10 m at the edge. The plots in figure 3.6 shows good agreement. From here it was concluded that quadratic formulation of the elements was most efficient in terms of computation time and precision in all cases. The chosen soil element types are henceforth C3D10, and CIN3D8 for the boundary region.

0.0 0.5 1.0 1.5 2.0

a0 1

2 3 4 5 6

Kev

Kobori et al. (1971) Present solution

0.0 0.5 1.0 1.5 2.0

a0 1

2 3 4 5

Cev

Figure 3.6: Comparison to solution by Kobori et al. (1971): Real and imaginary part of the impedance function.

3.3.2 Verification against a BEM model

To investigate the difference between using a structural damping and of viscoelastic damping, the model was compared to a BEM model with the same settings as for Kobori et al. (1971), but with a rigid plate and a loss factor η = 0.04 instead. Dimensionless frequency, a₀ = ωB/c_s, is used, and the real and imaginary parts are normalized with the static stiffness K0.

H E, ρ, ν, η H E, ρ, ν, η

F/m² F

Figure 3.7: Assumed conditions for verification according to BEM model with a rigid plate.

With the plot for the BEM model as the key, it is obvious that the difference between the use of viscoelastic damping and structural damping is small, see figure 3.8. Since structural damping is more convenient, it is used henceforth.

0.0 0.5 1.0 1.5 2.0

a0 0.2

0.4 0.6 0.8 1.0

Re(Z)/K0

0.0 0.5 1.0 1.5 2.0

a0 0.0

0.4 0.8 1.2 1.6

Im(Z)/K0

BEM

FEM Viscoelastic FEM Structural

Figure 3.8: Comparison to BEM solution: Real and imaginary part of the impedance function.

3.3.3 Verification against Zheng et al. (2014)

Zheng et al. (2014) [17] is a numerical solution of a 10 meter deep viscoelastic, homogeneous and isotropic soil with a vertical, one-dimensional elastic pile, with circular cross section, and is fixed to a rigid bottom. The pile is subjected to a harmonic vertical unit load and only impedance functions for vertical displacement is examined, see figure 3.9. Impedance functions are presented as functions of the non-dimensional frequency a = ωH/cp where cp is the S-wave velocity of the pile, and is dependent on its elastic modulus E_p and density ρ_p. The model according to Zheng et al. (2014) presents results for frequencies up to 500 Hz, but here simulations up to 200 Hz were performed. The impedance is rewritten in the non-dimensional form:

K_d= ZH

E_pA (3.2)

where A is the cross section area of the pile. The other parameters are:

r = 0.2 m H = 4 m E_p= 25 GPa E_s = 0.28 GPa ρ_p = 2500 kg/m³ ρ_s = 2200 kg/m³ η = 0.02 ν_p = 0.2 ν_s = 0.4

Es, ρs, νs, η

E_p, ρ_p, ν_p H

r F

Figure 3.9: Conditions for verification according to Zheng et al. (2014).

This gives a shear wave velocity of cs≈67 m/s. The structural damping was used so that ξeq= 1

%. The radius of the mid region was set to 2 m and had a mesh size of 0.1 m, while the total radius was 200 m with an assigned mesh size of 10 m at the edge. The piles were modeled both using embedded beams and 3D solids, with features discussed in section 2.2.4. Element type for the embedded beams were B31, C3D10 for the soil and CIN3D8 for the boundary region. Assigned element size for the pile was 0.1 m for both embedded beam and solid.

The plots in figure 3.10 are showing a comparison to the solution by Zheng et al. (2014) and between the two different modeling options for the piles.

0 1 2 3 4

a

−20

−10 0 10 20 30

Re(Kd)

0 1 2 3 4

a

0 10 20 30 40 50 60

Im(Kd)

Zheng et al. (2014) Embedded beam 3D solids

Figure 3.10: Comparison to solution by Zheng et al. (2014): Real and imaginary part of the impedance function.

Since agreement for the embedded beam option only seemed to be valid for frequencies up to a ≈2, it was of interest to study the lower frequency range more closely. This is shown in figure 3.11.

0.0 0.1 0.2 0.3 0.4 0.5

a

0.96 0.98 1.00 1.02 1.04

Re(Kd)

0.0 0.1 0.2 0.3 0.4 0.5

a

0.00 0.02 0.04 0.06 0.08 0.10 0.12

Im(Kd)

Zheng et al. (2014) Embedded beam 3D solids

Figure 3.11: Details of the lower frequency range in 3.10

As demonstrated, agreement was not seen for the embedded beam option even for very low fre- quencies. It appears more reliable to use solid elements for all frequencies.

Figures 3.12 and 3.13 are showing the converged three-dimensional solids solution with increased frequency resolution. Ripple effects are seen, probably due to the low damping in the soil (1

%). This indicates that further improvements of the model could have been made, primarily by increasing the mid-region radius. However, this was not further investigated.

0 1 2 3 4

a

−4

−2 0 2 4 6

Re(Kd)

0 1 2 3 4

a

0 2 4 6 8 10 12

Im(Kd)

Zheng et al. (2014) 3D solids

Figure 3.12: Comparison between FEM solution using 3D solids and solution by Zheng et al. (2014): Real and imaginary part of the impedance function, full range.

0.0 0.1 0.2 0.3 0.4 0.5

a

0.96 0.98 1.00 1.02 1.04

Re(Kd)

0.0 0.1 0.2 0.3 0.4 0.5

a

0.00 0.02 0.04 0.06 0.08 0.10 0.12

Im(Kd)

Zheng et al. (2014) 3D solids

Figure 3.13: Comparison between FEM solution using 3D solids and solution by Zheng et al. (2014): Real and imaginary part of the impedance function, lower frequency range.

More results of the convergence analysis for Zheng et al. (2014) is shown in Appendix A.

3.3.4 Verification against Padrón et al. (2012)

Similar to the Zheng et al. (2014) model, the Padrón et al. (2012) [18] model is a numerical solution with a pile group with either 2 × 2 or 3 × 3 inclined piles, with circular cross sections, in a homogeneous isotropic viscoelastic soil. The soil is solved by the boundary element method (BEM) and the piles are modeled using finite elements. The pile toes are fixed in bedrock, and the pile heads are fixed into a conceptualized stiff massless pile cap that can rotate, but is not in contact with the ground surface. A harmonic unit load is placed in the pile cap and the impedance function is registered in the origin of the ground surface, see figure 3.14. The unit load