3D dynamic response of submerged floating tunnels under seismic and hydrodynamic excitation

Engineering Structures 30 (2008) 268–281 www.elsevier.com/locate/engstruct

3D dynamic response of submerged floating tunnels under seismic and hydrodynamic excitation M. Di Pilato a , F. Perotti a,∗ , P. Fogazzi b a Department of Structural Engineering, Politecnico di Milano, Milan, Italy b Department of Civil Engineering, Architecture, Land Planning and Environment, Universit`a di Brescia, Brescia, Italy

Received 13 December 2005; received in revised form 29 March 2007; accepted 2 April 2007 Available online 8 May 2007

Abstract In the paper a numerical procedure is described for the dynamic analysis of seabed anchored floating structures, with particular reference to the so-called Archimedes bridge solution for deep water crossings; attention is devoted to the design solution encompassing slender bars as anchor elements. A geometrically nonlinear finite element, developed in previous work, is here refined extending its capabilities to full 3-D analysis and to nonlinear modelling of hydrodynamic loads due to steady current and wind waves. The element is implemented in a numerical procedure for the dynamic time domain step-by-step analysis of nonlinear discretized systems; consistently, hydrodynamic and seismic loading are introduced by generating artificial time-histories of spatially variable seismic motion and wind waves. An example of on application is shown regarding the behavior of the dynamic model of a submerged tunnel proposed for the Messina Strait crossing. The model is subjected to an extreme multiple-support seismic loading having a PGA equal to 0.64g and to an extreme wave loading with significant wave height of 16 m. The dynamic behaviour in the two loading situations is illustrated and compared, showing interesting facets, especially in terms of interaction between the tunnel and anchoring bars oscillations. c 2007 Elsevier Ltd. All rights reserved.

Keywords: Nonlinear dynamics; Geometrical effects; Hydrodynamic loads; Multiple-support seismic excitation

1. Introduction The submerged floating tunnel (SFT or Archimedes bridge) is an innovative concept for crossing waterways, such as sea-straits, fjords or lakes [1]. It basically consists of a tunnel floating in the water, due to positive net buoyancy, at some convenient depth and fixed to the seabed by means of a suitable anchoring system, encompassing cables or bars. Economical and environmental reasons support the belief that the Archimedes bridge can be an advantageous solution to the crossing problem in many circumstances. In spite of this, a first realization of a SFT is still missing, while a large number of more traditional immersed tunnels (i.e. tunnels directly supported by the seabed) are present worldwide. In terms of anticipated structural behaviour, the most striking difference between the two structural types seems to be the fact that the ∗ Corresponding author. Tel.: +39 223994229; fax: +39 223994220.

E-mail address: [email protected] (F. Perotti). c 2007 Elsevier Ltd. All rights reserved. 0141-0296/$ - see front matter doi:10.1016/j.engstruct.2007.04.001

SFT is much more prone to undergo significant oscillations due to dynamic effects (see [2,3]). These can be caused by different phenomena, such as: • supports motion due to seismic events; • water waves due to wind, to seismic motion of the seabed, or to water density effects (“internal” waves); • vortex shedding due to steady current; • dynamic forces due to traffic loading. Nevertheless the SFT solution have been proposed for crossings characterized by severe seismic and/or hydrodynamic design loads, such as the Messina Strait in Italy [4–6] the Funka Bay in Japan [7] and the Hogsfjord in Norway [8]. This justifies research efforts towards the development of numerical procedures for simulating the dynamic response of SFTs. Note that the need for such procedures is not exclusively related to exceptional loading: a favourable vibrational behaviour in normal operational conditions seems to be mandatory for the acceptance by the users of an infrastructure of this type. In fact, the vibration levels that are felt as acceptable on normal bridges

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

and footbridges would make people feeling insecure travelling, surrounded by water, through the most innovative construction of modern civil engineering. Though general purpose nonlinear structural analysis computer codes can be used, for example, for computing the seismic response of SFTs, the development of ad hoc numerical tools seems to be justified and advisable for covering the above listed dynamic loading conditions. This is motivated by numerical efficiency considerations and by the necessity of developing special loading models and analysis procedures; the cases of vortex induced vibration is typical example in this respect. In this light, in previous research work [9] the development of a numerical procedure for studying the nonlinear behaviour of Archimedes bridges under seismic and hydrodynamic loading was initiated; as a first step, attention was mainly devoted to the modelling of the anchoring system, whose behaviour is of extreme importance in determining the overall structural dynamic response. The design solution based on the use of slender bars for connecting the tunnel to the foundations was addressed and an ad hoc finite element was developed for efficiently modelling these elements; linearized hydrodynamic loading due to waves and current was introduced, within the framework of the Morison approach (see [10,11]). The analysis was restricted to oscillations in the transverse plane, i.e. in the plane which is normal to the tunnel axis. In the work described here the capabilities of the finite element and of the analysis procedure are extended to represent full 3D motion of the entire SFT. Modelling of the wind waves field has also been improved by introducing directional effects and hydrodynamic loading is now treated in the complete nonlinear form. The formulation of the finite element developed for anchor bars will be described in Section 2 of the paper, while Section 3 will be devoted to numerical procedures for response computation. In Section 4 an example of application will be described, regarding one of the proposed solutions [4] for the Messina Strait crossing. 2. The anchor element model (NWB element) A finite element to simulate the dynamic behaviour of an entire anchor element was developed; since for deep water crossings anchor bars can have extremely high slenderness values (up to 600 in the case here studied), particular attention was devoted to geometrical effects. A formulation based on the definition of a local moving reference, such as in the classical “corotational” approach [13], and on the introduction of relative coordinates to model the local transverse oscillations of the bar was devised. The formulation is essentially based on the hypotheses of hinged element ends and constant axial force. 2.1. Computation of the stiffness matrix The element stiffness matrix was computed under the hypothesis of linear elasticity and by assuming that transverse

269

Fig. 1. Anchor element: local coordinate system.

displacements and rotations relative to the bar chord are small; axial deformation, which is assumed to be small and constant along the bar, is expressed as a quadratic function of nodal displacements, including the effect of transverse motion. 2.1.1. Definition of local stiffness A local coordinate system is defined (Fig. 1) such that the i axis is attached to the element chord while the j and k directions are aligned to the principal axes of the element section. In the local system the bar ends are fixed, so that the element is statically determined. A simple model of its deformative behavior is obtained by introducing the following local coordinates: 1. v1 is the relative end displacement due to elastic effects, i.e. εi = 0i = v1 /L is the constant axial deformation; 2. v2 is the midspan transverse displacement in the j direction, i.e. v2 = v j (L/2); 3. v3 is the midspan transverse displacement in the k direction, i.e. v3 = vk (L/2). All displacements are assumed to be small; v1 leads to a small axial strain while v2 and v3 lead to curvature distributions which can be approximated by the second derivatives v 00j vk00 . Note that v1 , v2 and v3 can be seen as generalized small strain components for the element. By introducing sinusoidal shape functions the transverse displacement and curvatures can be expressed as  π 2 πξ πξ v2 → χ j ∼ sen v2 = v 00j = − L L L  π 2 πξ π ξ 00 vk (ξ ) = sen v3 → χk ∼ sen v3 . = vk = − L L L

v j (ξ ) = sen

(1)

By neglecting shear and torsional deformations the elastic stiffness matrix k can be defined, by standard manipulation, in the following way:   A 0 0 E 2  0 k =  0 J j π 4 /2L (2)  L 2 0 0 Jk π 4 /2L where A, J j and Jk are the area and the moments of inertia of the element section and E is the material Young’s modulus. The local stiffness k relates the coordinate vector v with local forces N (axial force) and F jt and Fkt (generalized forces associated to the midspan deflections), i.e.  T P = N F j Fk = kv. (3)

270

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

The introduction of nonlinear material behaviour can be performed, with standard criteria, at this stage of the formulation. In addition, the element kinematics can also be enriched by adding other deformation modes, for example, higher order harmonic shape functions. 2.1.2. Local to global coordinate transformation The global element coordinates u are defined such that (Fig. 2); u 1 , u 2 , u 3 , u 6 , u 7 , u 8 are the components of the total end displacements in the global reference system, while u 4 , u 5 are the midspan transverse displacements relative to the bar chord i.e. u 4 = v2 and u 5 = v3 . To relate the local coordinate v1 to the u displacements, the contributions due to ends displacement and to the effect (‘bowing effect’) of the transverse motion are treated separately and then superimposed. The first contribution can be expressed, retaining terms up to the second power, as " (e)

v1 = (u 6 − u 1 )α1 + (u 7 − u 2 )α2 + (u 8 − u 3 )α3 (u 1 − u 6 )2 (u 2 − u 7 )2 (u 3 − u 8 )2 + + + 2L 2L 2L

# (4)

α1 , α2 , α3 being the director cosines of the bar chord. By use of (1), the second-order effect of the relative transverse displacement is given by π2 2 = (u + u 25 ). (5) 4L 4 By adding the two contributions the total elongation can be (e) (b) obtained as v1 = v1 + v1 ; in incremental form the same quantity can be written as (b) v1

1 1v1 = v1 (u + 1u) − v1 (u) = 1u B L + 1uT B N L 1u (6) 2 where 1u is the vector of nodal displacement increments and where the linear and nonlinear generalized strain–displacement matrices are defined as follows: T  π2 π2 (7) BL = R S T u4 u 5 −R −S −T 2L 2L BN L 1 0 0 0 0 −1 0 0 T

0 0  1 0 =  L0  −1  0 0

1 0 0 0 0 −1 0

0 1 0 0 0 0 −1

0 0 π 2 /2 0 0 0 0

0 0 0 π 2 /2 0 0 0

0 0 0 0 1 0 0

−1 0 0 0 0 1 0

0 −1   0  (8) 0  0  0 1

u 3 −u 8 u 2 −u 7 6 being R = −α1 + u 1 −u L , S = −α2 + L , T = −α3 + L .

2.1.3. Computation of element stiffness matrix in the global reference system To compute the element tangent stiffness matrix to be used in a time integration procedure, the restoring forces at time t

Fig. 2. Anchor element: local to global coordinate system.

are supposed to be known. To estimate the same forces at time t + 1t the principle of virtual displacements can be written by assuming displacement increments 1u as well as axial strain and curvature increments (1εi , 1χ j , 1χk ) as kinematics quantities. By equating the internal work to the one done by generalized elastic forces Rt+1t , one can write, in the small deformation hypothesis Z δ (1u)T Rt+1t = δ (1ε N ) N t+1t dξ L Z
+ δ 1χ j M t+1t dξ j L Z + δ(1χk )Mkt+1t dξ. (9) L

Axial strain and curvature increments can be written, according to Eqs. (1) and (6), as   1v1 1 1 1ε N = = 1uT B L + 1uT B N L 1u L L 2 π2 π2 πξ πξ (10) = −1u 4 2 sin sin 2 L L L L π2 π2 πξ πξ = −1u 5 2 sin 1χk = −1v3 2 sin L L L L where the u coordinates appearing in the B L matrix are the values at time t. In Eq. (9), the contribution Wi,a of axial elongation to virtual work can be expressed in the following form: Z
1 Wi,a = δ (1v1 ) N t + 1N dξ L L Z Z 1 1 = δ(1v1 )1N dξ + δ (1v1 ) N t dξ. (11) L L L L

1χ j = −1v2

The first term appearing at the r.h.s. of Eq. (11) can be estimated by adopting the linear approximations δ(1v1 ) ∼ = δ(1uT )B L ;

1N = k11 1v1 ∼ = k11 BTL 1u. (12)

As to the second term, upon substitution of Eq. (6) we obtain:   Z Z  1 1 δ 1uT N t δ (1v1 ) N t dξ = B L dξ L L L L  Z 1  T t T + δ 1u N B N L dξ 1u . (13) 2 L

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

Finally, the flexural contribution to virtual work can be written, by use of Eqs. (1) and by incremental decomposition of the two bending moments, as Z δ(1χ j )(M tj + 1M j )dξ Wi, f 1 = L

Wi, f 2

= (k22 1u 4 ) δ (1u 4 ) + F jt δ (1u 4 ) Z
δ (1χk ) Mkt + 1Mk dξ =

(14)

L

= (k33 1u 5 ) δ (1u 5 ) + Fkt δ (1u 5 ) where F jt and Fkt are the generalized components of the restoring forces w.r.t. u 4 and u 5 . Upon substitution in (9) of Eqs. (12)–(14) the generalized restoring force at time t + 1t can be estimated through the equation Rt+1t = (K L + K N L ) 1u + Rt

(15)

where KL =

k11 L

KN L = Rt =

0

Nt L t Z

N L

L

Z

(16)

L

Z

0

B L BTL dξ + K L , f B N L dξ = N t B N L

0 L

 B L dξ + 0

0

0

(17) F jt

Fkt

0

0

0

T

(18)

are respectively the elastic, configuration dependent, stiffness matrix K L , the geometrical stiffness matrix K N L , while Rt lists generalized restoring forces at time t. Owing to the performed linearization, the restoring forces (15) are not in equilibrium with other loads. To enforce equilibrium a modified Newton–Raphson iteration is employed in the procedure here described (Section 3); this requires computation of actual elastic forces at time t + 1t. To do this, local coordinates can be evaluated at time t + 1t according to (4) and (5) and considering that u 4 = v2 and u 5 = v3 . Local internal forces can be subsequently computed by (3); once local forces are computed, generalized forces can be obtained upon substitution in Eq. (18). 2.2. Computation of the element mass matrix and static nodal loads To express the mass matrix the bar displacement field is linearized about its underformed configuration. The element absolute velocity field can be readily expressed in a new ¯ encompassing nodal displacement in coordinate system u, the local directions i, j, k (u¯ 1 , u¯ 2 , u¯ 3 , u¯ 6 , u¯ 7 , u¯ 8 ) and midspan relative displacements u¯ 4 = u 4 and u¯ 5 = u 5 , i.e.   ˙v¯ i (ξ ) = 1 − ξ u˙¯ 1 + ξ u˙¯ 6 L L   ξ ξ πξ v˙¯ j (ξ ) = 1 − u˙¯ 2 + u˙¯ 7 + sen u˙¯ 4 (19) L L L   ξ ˙ ξ πξ v˙¯ k (ξ ) = 1 − u¯ 3 + u˙¯ 8 + sen u˙¯ 5 . L L L

271

By writing the bar kinetic energy the local element mass matrix m can be defined as Z   1 1 T ˙ 2 2 2 ¯ (20) ρi v˙¯ i + ρ j v˙¯ j + ρk v˙¯ k dξ = u˙¯ mu. T = 2 L 2 In (20) the mass per unit length values, defined as the product of density times section area, have been introduced. Note that different values are used in the i, j, k directions to account for added mass effects. To get the element mass matrix M in the global reference frame a standard transformation accounting for the rotation of local reference system is subsequently performed on the coordinates representing the end displacements of the bar. Finally, consistent static nodal loads have been derived, under the same hypotheses adopted for inertia loads, to take account of element weight and buoyancy. The element formulation here described has proved to be very efficient: in [9] it was shown how to get the same results with a “general purpose” software package a finite-element discretization encompassing four elements per anchor bar was necessary and that the analysis procedure had to account for both large displacement and large strain conditions. Summing up, the proposed element was able to outperform, in terms of numerical efficiency, the standard procedure in the dynamic analysis of a tunnel section. The extension of the proposed finite-element capabilities to full 3-D motion is also of utmost importance since, as it will be shown in the following, both transverse and longitudinal oscillations (with respect to tunnel axis) of the anchor bars show critical aspects. 2.3. Element loads due to hydrodynamic effects The computation of hydrodynamic forces on a floating tunnel is one of the most difficult tasks in the simulation of the dynamic response under extreme and serviceability conditions, since it involves the complexity of the fluid–structure interaction phenomena [10]. Consistently with the adopted simplified approach to structural modelling, simple models were also sought for the hydrodynamic forces due to wind waves; these models were based on the Morison equation, whose application is largely justified for the anchor bars, given their diameter range. The extended approach given by Chakrabarti for the case of inclined elements [11] was used. Accordingly, the components along the global axes x, y, z of the wave force on a moving cylinder were written in terms of the normal components (w.r.t. the element axis) of relative velocity and acceleration, i.e. π D2 1 ρC D D|W |vx + ρC M w˙ x 2 4 2 πD − ρ(C M − 1) u¨ 4 1 π D2 = ρC D D|W |v y + ρC M w˙ y 2 4 π D2 − ρ(C M − 1) v¨ 4

fx =

fy

(21a)

(21b)

272

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

1 π D2 ρC D D|W |vz + ρC M w˙ z 2 4 π D2 w. ¨ − ρ(C M − 1) 4 In Eqs. (21) the following quantities are introduced:

fz =

(21c)

ρ, D, C D , C M are the water density and the element diameter, drag and inertia coefficients: u, v, w are the projections along x, y, z of the normal component of the element displacements; wx , w y , wz are the projections along x, y, z of the normal component of the wave velocity obtained as the difference between the total velocity components u x , u y and u z and the tangential velocity components, i.e. w x = u x − α1 u i ; w z = u z − α3 u i

w y = u y − α2 u i ;

where α1 , α2 , α3 are the director cosines of the bar and u i is the tangential velocity: u i = u x α1 + u y α2 + u z α3 vx , v y , vz are the projections of the fluid–structure relative velocity, expressed as vx = wx + Ux − u; ˙ vz = wz + Uz − w˙

v y = w y + U y − v; ˙

where Ux , U y , Uz are the components of the steady current velocity. |W | is  the modulus of the relative velocity, i.e. of the vx v y vz vector. The first term on the right-end side of Eq. (21) represents the drag loading, the second the inertia loading, the third the added mass effect. The distributed load (21) has been treated by means of a lumped parameter approach: to do this transverse concentrated loads acting at generic point P, representing either the element ends (point A and C) or the midspan (point B), have been computed as follows: fP

=

1  ρC D A P |W P | v P x + ρC M V P w˙ P x − ρ (C M − 1) V P u¨ P   2

    

D P P Py M P Py M P P 2     1 ρC A |W | v + ρC V w˙ − ρ (C − 1) V w¨ D P P Pz M P Pz M P P

   

1

ρC A |W | v

2

+ ρC V w˙

− ρ (C

− 1) V v¨

(22)

where A p and V p are the “relevant” area and volume at element nodes A, B or C (P = A, B or C). For a circular section and introducing a parameter ϕ, that defines the area and volume portions which are relevant for the midspan point, we can write DL ; A B = ϕ DL; 2 π D2 π D2 V A = VC = (1 − ϕ) L; VB = ϕ L. 8 4 The generalized components of the forces (22) referred to the element co-ordinates u are the equivalent nodal A A = AC = (1 − ϕ)

dynamic loading and added mass coefficients. They have been subsequently computed by means of standard manipulation based upon virtual work done by the force themselves. The described model of wave actions leads, on one hand, to nonlinear dynamic forces on the structure and, on the other hand, to a modification of the mass matrix by terms that take account of the water motion. An identical procedure was set for defining the hydrodynamic forces on the tunnel, even though, for typical diameters of proposed SFTs, the Morison approach can be justified only for extremely large waves. For waves having larger probability of occurrence, which are of interest for fatigue and serviceability considerations, diffraction effects must be considered: in this case, however, given the lower excitation intensity, linear analysis is justified on structural grounds. Diffraction forces on a moving body, on the other hand, admits a linearized solution encompassing frequency-dependent added mass and hydrodynamic damping effects (see [10]) which can be efficiently handled, in a totally linear setting, by direct frequency domain analysis. Such an approach is under investigation by this research group, for studying the tunnel response under “normal” waves and traffic loading. 3. Equations of motion and solution strategies The equations of motion of a soil–structure system subjected to multiple-support seismic excitation as well as to other dynamic loads (e.g. hydrodynamic effects) can be written in the matrix form Mq¨ + Cq˙ = R + Qs + Qh + Q

(23)

where q is the vector of Lagrangian coordinates (here representing total generalized displacements), M and C are the inertia and damping matrices, while R, Qs , Qh and Q are vectors listing, respectively, the generalized components of the nonlinear restoring forces, of the “equivalent” seismic forces, of the nonlinear hydrodynamic forces (22) and of other static or dynamic actions. At time t the restoring force vector R sums generalized components of both elastic forces and elemental nonlinear forces (18). If we assume linear behaviour of the ground and lumped-parameter (frequency independent) modelling of soil–structure interaction, the seismic forcing term can be expressed in the following partitioned form:     0 0 Qs = + (24) (g) ( f ) (g) ( f ) Ccc q˙ c Kcc qc (f)

(f)

where q˙ c and qc are the vectors of free-field ground velocity and displacement at soil–structure “contact” points, i.e. at the (g) points located at the soil–foundation interface, while Ccc and (g) Kcc are respectively the soil damping and stiffness matrices referred to “contact” degrees-of-freedom as well. For the formation of the damping matrix accounting for structural dissipation, the procedure described in [9] was adopted: accordingly, a damping matrix is first formed for a reduced structural model of the SFT in which anchor bars are

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

represented by a single truss element, thus disregarding local transverse oscillations. To build this matrix, the iterative method proposed in [9] was used: the method allows for the formation of a banded matrix, though preserving a good approximation of the damping factor for a significant number of normal modes of the linearized model. The matrix is subsequently expanded by adding rows and columns corresponding to the coordinates accounting for local oscillations of the anchor elements. Given the definition of these coordinates, which are relative to the motion of the element chord, damping of the oscillations can be directly obtained by adding suitable coefficients on the corresponding elements of the damping matrix diagonal. For each component of transverse motion (local directions q j and k) the damping coefficient is

simply computed as 2ς K ll0 Mll where l is the DOF number in

the global structural model, ς is the damping factor and K ll0 is the initial value of the transverse local stiffness. Note that ς can be different from the damping factor adopted for the “global” tunnel oscillations. It must be observed that the described procedure for efficient and consistent definition of the damping matrix is made possible by the introduction, here proposed and developed, of relative coordinates for describing local transverse oscillations of the anchor bars, which again demonstrates to be highly advantageous for the problem at hand. Note that a similar situation arises in apparently different structural systems, such as large cable-stayed bridges. Finally, the Eq. (23) were solved via a numerical procedure based upon the Newmark method for time integration and the Modified Newton Raphson (MNR) scheme for equilibrium iteration at each time step. The procedure [14,17] is suitable for treating both nonlinear restoring forces and applied dynamic forces showing nonlinear dependence on the system velocity and configuration, the latter being the case of drag forces (22). Response estimation within a time step is performed, in the Newmark procedure, upon linearization of the restoring force increments as expressed, at the element level, by (15)–(17). The restoring forces contribution to the iteration matrix for the MNR scheme is obtained as well by assembling elemental matrices (16) and (17). On the other hand, the adopted iterative strategy [14] does not encompass the formation of fluid-dynamic damping and stiffness matrices. Extensive numerical testing has demonstrated that the resulting slight increment in the number of iterations at each timestep, which lies anyway in the range of some units, is largely compensated by the advantage of avoiding the treatment of nonsymmetrical matrices. Note, in addition, that the effect of configuration change on hydrodynamic loading, to which the use of a fluid-dynamic “load stiffness” matrix is related, is here disregarded since the normal component of fluid velocity appearing in expressions (21) is computed with respect to the initial position of the bar. For a more refined analysis, which is necessary for more flexible elements, the reader is referred to the criteria given in [14], which are specific to the case of cables, or to the more general and comprehensive treatment offered in [15].

273

4. Environmental actions modelling and simulation The environmental action fields, i.e. the seismic motion field and the sea wave–current field were both introduced by simulation techniques: in both cases artificial time-histories were generated by the wave superposition technique, according to the stochastic models presented in [10] (sea waves) and [21, 22] (seismic motion). The adopted models for hydrodynamic loading are consistent to the hypothesis of a crossing in deep water, such as the one in the example here shown. Within this context, simple models have been chosen for the wave kinematics and for the probabilistic modelling of the wave field. More sophisticated assumptions, however, can be easily introduced into the artificial simulation process. 4.1. Sea waves and current A random ocean wave model, the Multivariable Stationary Model [10], was adopted within a simulation technique based on the superposition of waves having the form η(x, z, t) =

N X M X

αmn cos[(km cos φn )x

n=1 m=1

+ (km sen φn )z − ωm t + θmn ].

(25)

The only random variable is here the phase angle θmn , having uniform probability density function within the interval (0, 2π), while the φn angle defines the wave direction and the parameter km refers to the wave number. The circular frequency ωm is in turn related to km by the dispersion condition. The amplitudes αmn were determined introducing the so-called directional wave spectral density S (ωm , φn ) [10] as a function of two variables, the frequency ωm and the wave direction φn . This spectral density function is defined as S (ω, φ) = S (ω) G (ω, φ)

(26)

where S(ω) is the power spectrum of the wave height at a point and G (ω, φ) is a directional spreading function. In the first test applications of the procedure a very simple, frequencyindependent form of the directional spreading function has been introduced, namely the following “cosine-squared spreading function”: (2 π cos2 φ |φ| < (27) G(φ) = π 2 0 elsewhere. In a more refined analysis the Mitsuyasu [12] spreading function, accounting for frequency dependence of the directional effect, can be easily implemented in the artificial generation procedure. The power spectrum here adopted is the Pierson–Moskowitz one, having the well-known form "  −4 # ω 2 −5 (28) S (ω) = αg ω exp −1.25 ω0

274

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

where α = 0.0081 and the ω0 frequency is related to the significant wave height Hs through the relation ω02 = 0.161g/Hs . At each frequency in the summation (25) the variance, i.e. the area under the spectral curve, provides the amplitude of the wave elevation according to the relation 1 2 α = S ( f m , φn ) 1 f 1φ. (29) 2 mn To examine the influence of currents on a random waves field, a very simple modification to directional wave spectral density was introduced according to [16], in the following equation:   1 S ( f, φ) = S ( f ) G ( f, φ) · (30) 1 + 2ωr |U |/g where ωr = ω − K T U. U and K are, respectively, the current velocity vector (U = Ux Uz ) and wave number  vector (K = k cos φ k sin φ ). The modified spectrum (30) expresses the physical occurrence that when a wave encounters a current, wave characteristics undergo changes. In a following current, where the flow of the current is in the direction of wave propagation, the wave amplitude decreases, and the wavelength increases so that the sea becomes calmer. On the other hand, in an adverse current, where the current opposes the wave, wave amplitude increases and wavelength decreases. For each component in (25), velocity and acceleration frequency components were subsequently expressed according to Airy linear wave theory [11]. Note that more sophisticated theories can be rather easily handled by the artificial simulation procedure; given the water depth, however, linear theory is justified in the example here shown. 4.2. Seismic motion The integration of Eq. (23) requires the knowledge of seismic motion time histories, in free-field conditions, at all soil–structure interface nodes. To this aim a numerical procedure for artificial ground motion generation has been adopted [18], which is based upon the spectral representation method of Shinozuka [19] and is developed in accordance with the criteria reported by Monti et al. in [20]. Simulation is based on a stochastic model; the free-field surface acceleration at a point a(t) is modelled as an uniformly modulated non stationary process. This can be expressed as the product of an ‘underlying’ stationary process a¯ (t) and of a deterministic envelope η (t). To describe the underlying stationary acceleration at a site, the Kanai-Tajimi power spectral density (PSD) as modified by Ruiz and Penzien [21] is adopted, expressed as
2 1 + 4ςg2 ω/ωg Sa¯ (ω) = S0 · h
2 i2
2 1 − ω/ωg + 4ςg2 ω/ωg
4 ω/ω f ·h
2 i2
2 1 − ω/ω f + 4ς 2f ω/ω f

(31)

where S0 is a scale factor, ωg and ςg are the characteristic ground natural frequency and damping, while ω f and ς f are the parameters of a high pass filter introduced to avoid singularities, in the lower frequency range, of velocity and displacement amplitudes. The same model is assumed for all contact points, which means that global attenuation and site effects are disregarded. Under such hypothesis, to describe the spatial variation of the ground motion the cross-spectral density between stationary acceleration at two sites (say P and Q) located at a ξ separation distance, have been expressed as the product Sa¯ P Q (ω, ξ ) = γ P Q (ω, ξ ) Sa¯ (ω)

(32)

of the PSD and of the coherency function γ P Q (ω, ξ ); for the latter a model was adopted which is able to represent both the ‘geometrical incoherence’ effect and the ‘wave passage’ effect. The first one is caused by reflections and refractions of waves along the propagation path, while the second, related to propagation velocity, affects the arrival time of wave trains at different sites. The chosen model is due to Luco and Wong [22] and is defined by the following expression: γ P Q (ω, ξ ) = exp[− (ωαξ/vs )2 ] exp[iωξ (L) /vapp ]

(33)

where vs is the shear wave velocity, vapp is the surface apparent velocity of ‘dominant’ waves, ξ (L) is the projection of the horizontal distance between the sites along the propagation direction and α is a parameter controlling coherence decay. Note that the apparent wave velocity should be, in principle, frequency dependent. In the model here adopted, however, constant apparent propagation velocity is assumed, so that the wave passage effect (complex factor in expression (33)) simply results in a delay of the arrival time of the ground shaking at the different structural support points. Since this can be easily handled by the structural analysis procedure, the real factor only, related to geometrical incoherence, has been retained in Eq. (33) for the artificial generation. Once a set of stationary acceleration segments has been generated according to expressions (31)–(33), to obtain the complete non-stationary ground motion field the following steps are performed: 1. Scaling of each acceleration time history by the following deterministic envelope:  2 t η (t) = for 0 6 t 6 t1 t1 η (t) = 1 for t1 6 t 6 t2   t − t2 η (t) = exp ln β for t2 6 t 6 tmax . tmax − t2

(34)

2. Double integration to get time histories of ground velocity and displacement, performed according to the frequency domain procedure described in [23].

275

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

Fig. 3. Side and front view of the floating tunnel example.

5. Example of application 5.1. Structural and external actions modelling The capabilities of the procedure to simulate the SFT response in a real case were tested by studying the dynamic behaviour of a tunnel model (Fig. 3), which is largely based on published data and design criteria [4] about the crossing of the Messina Strait between Punta S. Ranieri and Catona (4680 m). The scope of the numerical experiment is to show how the modelling approach here proposed, though simplified in many senses, is capable to deliver an overall picture of the main aspects and problems related to the SFT dynamic behaviour under seismic and wave excitation. The model is based upon the following hypotheses and criteria: • The tunnel is modelled by means of 3-D beam elements, neglecting axial deformation and torsional behavior. The structural section is a composite steel–concrete thick pipe (external diameter 15.95 m, internal 13.95 m) having an equivalent area, referred to concrete elastic properties, of 58.24 m2 and an equivalent moment of inertia of 16 370 m4 . “Equivalent” means here that in computing section parameters the steel areas have been multiplied by the steel-to-concrete elastic modulus ratio. Lumped masses are positioned at each anchoring section, having a 72 m constant spacing. A weight per unit length of 1200 kN/m is considered for the tunnel, resulting in a net positive buoyancy (difference between buoyancy and weight) equal to 455 kN/m. • The anchoring system lies in the plane normal to the tunnel axis (x y plane see again Fig. 3) and is here modelled through a couple of equivalent NWB (Fig. 4(a)) elements having the same total area and slenderness of the real bars. Properties of the anchor bars are indicated in Table 1. • Soil–structure interaction effects are accounted for by inserting three elastic springs at the bottom of each anchor element together with linear dashpots, representing wave radiation effect, acting in parallel to springs. Stiffness and damping constants computation are based upon the hypothesis of a six-pile foundation block. Single piles properties have been first computed, by using half-space theory (see Table 2); group effect has been subsequently introduced by use of static interaction coefficients [24]. • Given that no information was available in [4] on the tunnel-shore connection, the same spring–dashpot elements as computed for anchor block foundations are inserted, as a first rough assumption, at each tunnel end.

Fig. 4. (a) Simplified model of anchor bar system; (b) Displacements components for a typical section modelled by the NWB approach. Table 1 Geometrical and mechanical characteristics of anchoring system: (D, s) diameter and thickness of the real element, (A, J ) area and moment of inertia of the equivalent anchor element Type

D (m)

s (m)

A (m2 )

J (m4 )

A B C

1.850 1.950 2.058

0.062 0.065 0.068

0.696 0.770 0.850

0.278 0.342 0.421

• A damping factor of 0.03 was assumed for the structural dissipation effects. Table 3 lists the natural periods of the first 100 normal modes of the reduced model, i.e. of the model (see Section 3) which does not account for transverse oscillations of the anchor bars. Coming to the model adopted for seismic loading, the parameters in Table 4 were supplied for the definition of the

276

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

Table 2 Foundation piles properties: (L/r ) pile slenderness; (G s /E p ) soil-to-pile stiffness ratio; µs soil Poisson coefficient; A, J pile area and moment of inertia; k and c stiffness and damping constants G s /E p

L/r

µs

40 250 0.25 Stiffness and damping constants Single Group

Vertical motion Horizontal motion Vertical motion Horizontal motion

A (m2 )

J (m4 )

4.91

1.92

k (kN m−1 ) 6 704 860 1 035 725 17 200 000 2 870 000

c (t s−1 ) 44 462 9614 114 000 26 600

acceleration PSD, of the coherency function and of the envelope function equations (32)–(34). The assumed S0 value is consistent with a peak ground acceleration of 0.64g; apparent wave propagation goes from Calabria to Sicilia. The same stochastic properties are assumed for horizontal and vertical seismic motion. Independent generation of three sets of time histories, one set for the vertical component and the others for the longitudinal and transverse ones (w.r.t. to tunnel axis), was performed. 5.2. Results of the analyses Some of the results of the performed analyses are shown in the following sections. For selected tunnel nodes (see Fig. 3), three displacement components (transverse x, vertical y and longitudinal z) are reported. For the anchor elements transversal displacements will be shown; we will assume

that local coordinate v2 lies in the xy plane and we will call it TT (transversal–transversal) local displacement, while v3 is parallel to the tunnel axis and will be termed as TL (transversal–longitudinal) component (Fig. 4(b)). The results are reported at three tunnel sections (see Fig. 3): extreme (sec. 67 or 1), intermediate (sec. 50 or 18) and midspan (sec. 34) sections. The results shown, though incomplete, seem to be sufficient for illustrating some aspects which are typical of the overall dynamic behaviour of the structural system and, at the local level, of the response of the anchor bars. For a more clear understanding of the effects of each loading condition, seismic motion and hydrodynamic loading were separately applied to the SFT model. 5.2.1. Seismic excitation In the case of the seismic motion, three different loading conditions have been considered, in which the vertical (y), transverse (x) and longitudinal (z) components of seismic motion at a single point are supposed to be mutually uncorrelated. Two of these conditions correspond to partial correlation (expressions 32 and 33) of the ground motion at different points, the first one having a PGA of 0.64g in the transverse (x) direction and a PGA equal to 85% of this value applied in the longitudinal direction (z). In the second one the highest excitation acts in the longitudinal (z) direction (and 85% of the full excitation is applied along x while PGA in the vertical direction is always set to 0.64g. In the third seismic condition a synchronized and homogeneous ground motion is considered for comparison, where the same displacement and velocity time histories are imposed at all support points.

Table 3 Natural periods (s) of the reduced dynamic model (mode numbers and directions of motion are in italic) 1L 109.9 11 V 1.931 21 V 1.749 31 L 1.580 41 V 1.402 51 1.201 61 1.050 71 .9051 81 .7732 91 .6560

2L 4.710 12 T 1.899 22 T 1.712 32 T 1.552 42 TV 1.367 52 1.187 62 1.027 72 .8816 82 .7504 92 .6546

3L 2.365 13 V 1.899 23 V 1.712 33 V 1.552 43 TV 1.366 53 1.186 63 1.020 73 .8779 83 .7480 93 .6357

4T 1.970 14 T 1.869 24 TV 1.685 34 T 1.513 44 TV 1.332 54 1.158 64 .9942 74 .8522 84 Z .7480 94 .6331

5V 1.970 15 V 1.869 25 T 1.670 35 V 1.513 45 TV 1.328 55 1.150 65 .9902 75 .8512 85 .7243 95 .6156

6T 1.966 16 T 1.838 26 V 1.670 36 TV 1.475 46 TV 1.296 56 1.118 66 .9615 76 .8249 86 .7234 96 .6115

7V 1.966 17 V 1.838 27 T 1.630 37 TV 1.475 47 TV 1.286 57 1.114 67 .9608 77 .8249 87 .6997 97 .5964

8T 1.954 18 T 1.792 28 V 1.630 38 TV 1.438 48 TV 1.261 58 1.081 68 .9333 78 .8011 88 .6996 98 .5892

9V 1.954 19 V 1.792 29 T 1.590 39 TV 1.438 49 TV 1.244 59 1.081 69 .9324 79 .7989 89 .6770 99 .5781

10 T 1.931 20 T 1.749 30 V 1.590 40 T 1.403 50 TV 1.224 60 1.053 70 .9095 80 .7766 90 .6766 100 .5650

Table 4 Seismic motion modelling parameters ωg (rad s−1 )

ςg

ω f (rad s−1 )

ςf

α

vs (km s−1 )

vapp (km s−1 )

β

t1 (s)

t2 (s)

tmax (s)

10 s

0.4

1.0

0.6

1.0

3.5

4.95

0.25

2.0

10

40

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

277

Fig. 5. Displacement time histories (m) under seismic excitation: Tunnel at section 34.

Fig. 6. Displacement time histories (m) under seismic excitation: Tunnel at section 67.

Maximum transverse (x) and vertical (y) displacements of the tunnel are very similar for the three load cases (Figs. 5 and 6). The overall extreme is obtained for the case of maximum PGA in transverse direction and is close to 1 m; time-histories of these components show significant oscillations at the lower natural frequencies of the tunnel reduced model (see modes 4–11 in Table 3). These oscillations are obviously smaller at the tunnel ends (Fig. 6), where the time variation of the total displacement is dominated by ground motion. Longitudinal tunnel displacement time-histories reach their extreme values for the simultaneous excitation case and are characterized by lower frequency content when compared to transversal and vertical motions. Decreasing amplitudes occur getting closer to the end sections. The dominating oscillations occur at frequencies close to the ones of modes two and three of the linear model, while the first longitudinal mode, which involves an almost rigid body motion, is practically not excited.

In Fig. 7 time histories of relative transverse displacement components (TT or v2 and TL or v3 ) are shown. The superposition of low and high frequency oscillations can be easily noted in the behaviour of very slender elements (at section 34, Fig. 7(a)): the former are due to local transverse modes, while the latter occur at frequencies close to the ones of the relevant tunnel motion, which are very close to the natural frequencies of the reduced model. For longitudinal motion, resonance-type phenomena coupling local and global oscillations give rise to larger amplitudes at anchor bars having local transverse natural frequencies close to the ones of longitudinal tunnel modes. In Fig. 7(b) the phenomenon is shown for anchor bars at section 18, whose longitudinal (TL) flexural motion is close to resonant conditions with the second tunnel longitudinal mode (see Table 3). Short elements (e.g. section 66, Fig. 7(c)) show highfrequency small-amplitude transverse motion: this combines

278

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

Fig. 7. Relative displacement t.h. (m) under seismic excitation: (a) right anchor bar at sections 34; (b) left anchor bar at sections 18; (c) right anchor bar at sections 66.

with high axial forces, resulting in normal stress values which largely exceed the yield limit of normal constructional steel. A refinement of the procedure for including nonlinear material behaviour is therefore advisable. It can be also considered that modern seismic design is based on the reduction of acceleration response through controlled inelastic deformation in selected elements. In this respect the behaviour of short anchor bars, if properly designed and detailed, could be exploited to create “dissipative” mechanisms within the SFT. To get an overall picture of the above described phenomena, in Fig. 8 envelopes of extreme displacement response are given for non-synchronized seismic input having maximum intensity in the transverse (x) direction: for each section the figure gives peak values of tunnel displacement components and of local transverse (TT and TL) components of anchor

bars displacements. The curves clearly show the intensity of the “end-restraint” effect, giving rise to very high response of short anchor bars. The effect of resonance-type phenomena at sections from 10 to 20 is also clearly evidenced from longitudinal (TL) anchor bar response. On the other hand, extremes of tunnel displacement components show a very regular behaviour. Finally, when maximum excitation is applied in the longitudinal direction, the corresponding relative end displacements (tunnel/coast) reach a value of 0.72 m for Calabria coast and of 0.45 m for Sicilia coast. 5.2.2. Hydrodynamic loading due to waves and current The results presented in this paragraph are referred to the environmental conditions defined in [6] on the basis of previous

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

279

plane containing the tunnel axis) and steady hydrodynamic effects due to current (antisymmetrical). As can be inferred by initial (t = 0) TT displacement values, the two effects sum up for the right bar and subtract for the left one. Peak-to-peak extremes of relative transverse TT displacements reach amplitude values of about 2 m, close to those obtained in the seismic condition. Due to 3-D modelling of the wave field and to the lack of symmetry in the (z, y) plane, longitudinal displacements are significant (peak-to-peak exceeding 1 m), even though forces act prevalent in the transverse direction. 6. Conclusions

Fig. 8. Envelopes of displacement extremes under seismic excitation (a) TT components for anchor bars, transverse (x) and vertical (y) components for tunnel; (b) TL components for anchor bars and longitudinal (z) component for tunnel.

studies regarding the Messina Strait. The following data are common to all cases: • depth of tunnel relative to the mean wave profile equal to 40 m; • current velocity, acting in transverse direction (x), uniform with depth and equal to 1 m/s; • wave field having prevalent direction in the transverse direction. Extreme wave conditions are stated in [6] in terms of a significant wave height Hs = 13.5 m, which corresponds in the area to a return period equal to 100 years. Here results are shown for two values, namely Hs = 5 m and Hs = 16 m. In Fig. 9 the time histories of tunnel displacement components at a typical section are shown. It can be noted how response is characterized by significant harmonic components at frequencies which are lower than those obtained for seismic motion and that time histories reflect the aspect of velocity and acceleration of the wave/current field. Amplitudes obtained for the vertical and transverse displacement components of the tunnel are very similar, while the longitudinal component is obviously much smaller: displacement components’ extreme values are roughly one order of magnitude smaller than those obtained for the seismic load case. Components of relative transverse displacement (TT and TL) of the anchor rods of a typical section, illustrated in Fig. 10, show a nonsymmetrical behaviour of left and right elements, displacement amplitudes on the right side generally being larger. This behaviour is mainly influenced by the combination of buoyancy effects (symmetrical wrt, the vertical

A procedure for the nonlinear dynamic analysis of submerged structures under seismic or hydrodynamic excitation is being developed, especially to the aim of providing a reliable numerical tool for the design of innovative structures such as submerged floating tunnels. By the point of view of modelling, the key point of the procedure lies in the representation of the behaviour of the anchor elements connecting the tunnel to the seabed. The design choice here considered makes use of slender bars (of tubular section in the example here shown) as anchoring elements. The use of cables is under investigation in parallel research activity. Seismic analysis have been performed under the hypothesis of complete 3-D multiple-support excitation while response to hydrodynamic effects have been computed by considering nonlinear drag forces due to steady current and wind waves. An example has been analysed regarding a design proposal for the Messina Strait crossing, characterized by a total length of 4680 m and a maximum depth of about 285 m; extreme seismic loading was characterized by a PGA of 0.64g, while extreme hydrodynamic loading encompasses a significant wave height of 16 m and a steady current velocity of 1 m/s. The results here obtained point out some of the critical aspects of the tunnel behaviour under earthquake loading. When transverse oscillations are examined, tunnel behaviour appears to be rather predictable in terms of its linear behaviour. Flexural oscillations of long anchoring elements occur at much lower frequencies, while short ones, located at the tunnel sides, are subject to very high stressing, leading to inelastic behaviour. Longitudinal motion under seismic loading shows internal coupling phenomena between longitudinal tunnel modes and flexural oscillations of some anchor bars. This behaviour can be strongly affected by the tunnel end-restraints and by a possible different layout of the anchoring system; nevertheless the phenomenon seems to deserve further investigation and a particular care by a design standpoint. In fact, this internally coupled behaviour could be driven also by other excitation mechanisms, not considered here, such as vortex shedding due to transverse steady current. It can be noted, in this respect, that the low natural frequencies of transverse oscillations of anchor elements can be in the range of Strohual values, i.e. of frequencies governing the alternate shedding of vortices from the bar/cable section. A technique for simulation of

280

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

Fig. 9. Displacement time histories (m) under hydrodynamic excitation: tunnel at section 50.

Fig. 10. Relative displacement t.h. (m) under hydrodynamic excitation; right (top) and left (bottom) anchor bars at section 34.

vortex induced vibrations of the anchor bars is presently under investigation by the authors. Finally, given the tunnel depth and eigenspectrum, its dynamic behaviour under hydrodynamic loading due to

wind waves is characterized by very small amplitudes. Low frequency components of the wave field, however, are capable of exciting transverse flexural oscillations of the same order of magnitude as detected in extreme seismic conditions.

M. Di Pilato et al. / Engineering Structures 30 (2008) 268–281

Acknowledgments The work done by Sara Garavaglia and Simone Ferrario as a part of their graduation thesis is acknowledged. This research has been partially supported by MIUR (Ministry of Education, University and Research) under the project “Dynamic behaviour of structures: theory and experiments” (COFIN 01-02 & 03-04, website http://www.disg.uniroma1.it/fendis) The financial support of the “Ponte di Archimede S.p.A.” company is also gratefully recognized.

[10]

References

[15]

[1] FEHRL report no 1996/2a. Analysis of the submerged floating tunnel concept. In: FEHRL (Forum of European National Highway Research Laboratories). Berkshire (Crowthorne): Transport Research Laboratory; 1996. [2] Grantz WC. IT-SFT what’s the difference? Not really very much!. In: Proc. of the seminar on submerged floating tunnels held in Varenna. 2000. [3] Brancaleoni F, Castellani A, D’Asdia P. The response of submerged tunnels to their environment. Eng Struct 1989;11:47–56. [4] Bruschi R, Giardinieri V, Marazza R, Merletti T. Submerged buoyant anchored tunnels: Technical solution for the fixed link across the strait of Messina Strait crossings. In: Strait crossings. Rotterdam: Balkema; 1990. [5] Faggiano B, Landolfo Re, Mazzolani FM. Design and modelling aspects concerning the submerged floating tunnels: An application to the Messina Strait crossing. In: Proc. of the third international conference on strait crossings. 2001. [6] R.I.N.A. report. Verification of the basic design called Archimedes Bridge for the company Ponte di Archimede nello Stretto di Messina S.p.A. 1984. [7] Fujii T. Submerged floating tunnel in funka bay design and execution. In: Proc. of the international conference on submerged floating tunnels. 1996. [8] Skorpa L. Innovative norwegian fjord crossing. How to cross the Høgsjord, alternative methods. In: Proc. of the 2 ◦ congress A.I.O.M. (Marine and Offshore Engineering Association). 1989. [9] Fogazzi P, Perotti F. The dynamic response of seabed anchored floating

[11] [12] [13]

[14]

[16] [17]

[18]

[19] [20] [21]

[22] [23]

[24]

281

tunnels, under seismic excitation. Earthquake Eng and Struct Dyn 2000; 29:273–95. Sarpkaya T, Isaacson M. Mechanics of wave forces on offshore structures. Van Nostrand Reinhold Company; 1981. Chakrabarti SK. Hydrodynamics of offshore structures. Computational Mechanics Publications–Springer-Verlag; 1987. Mitsuyasu H. Observation of the directional spectrum of ocean waves using a cloverleaf buoy. J Phys Oceanogr 1975;5:750–60. Meek JL, Tan HS. Geometrically nonlinear analysis of space frames by an incremental iterative technique. Comput Methods Appl Mech Engrg 1984;47:261–82. Martinelli L, Perotti F. Numerical analysis of the non-linear dynamic behaviour of suspended cables under turbulent wind excitation. Int J Struct Stab Dyn 2001;1:207–34. Schweizerhof K, Ramm E. Displacement dependent pressure loads in nonlinear finite element analyses. Comput Struct 1984;18(6):1099–114. Hedges TS, Anastasiou K, Gabriel D. Interaction of random waves and current. J Waterway, Port, Coastal Ocean Eng (ASCE) 1985;111:275–88. Perotti F, De Amici A, Venturini P. Numerical analysis and design implications of the seismic behaviour of one-storey steel bracing systems. Eng Struct 1996;18(2):162–78. Martinelli L. Implementazione di un metodo numerico per la generazione di storie di accelerazione per strutture a grandi dimensioni in pianta. Ph.D. thesis (in partial fullfillment). Milano (Italy): Dipartimento di Ingegneria Strutturale, Politecnico di Milano; 1996 [in Italian]. Shinozuka M. Digital simulation of random processes and its applications. J Sound Vibration 1972;25(1):111–28. Monti G, Nuti C, Pinto PE. Non-linear response of bridges under multisupport excitation. ASCE J Struct Eng 1996;122(10):1147–59. Ruiz P, Penzien J. Probabilistic study of the behavior of structures during earthquakes. Berkeley (Calif.): Earthquake Engineering Research Center, University of California; 1969. UCBEERC-69/03. Luco JE, Wong HL. Response of a rigid foundation to a spatially random ground motion. Earthquake Eng Struct Dyn 1986;14:891–908. Perotti F. Frequency domain analysis of strong-motion accelerograms. Dipartimento di Ingegneria Strutturale del Politecnico di Milano, T.R. N. 3/87; 1987. Moore PJ. Analysis and design of foundations for vibrations. Rotterdam: AA Balkema; 1985.