structure interaction of submerged floating tunnels

Computers and Structures 72 (1999) 659±685

www.elsevier.com/locate/compstruc

Dynamic response and ¯uid/structure interaction of submerged ¯oating tunnels Svein Remseth*, Bernt J. Leira, Knut M. Okstad, Kjell M. Mathisen, Terje HaukaÊs Department of Structural Engineering, Norwegian University of Science and Technology, Trondheim, Norway Received 16 October 1998

Abstract Alternative approaches to stochastic dynamic response analysis of submerged ¯oating tunnels subjected to wave loading are presented. For the purpose of establishing force, damping and mass coecients for structural elements with three-dimensional ¯ow conditions, ¯uid/structure interaction is modeled as ®nite element implementation of the Navier±Stokes equation. The numerical examples emphasize the e€ects of wave direction, shortcrestedness, damping, geometrical sti€ness and frequency dependence in mass and damping coecients. # 1999 Elsevier Science Ltd. All rights reserved.

1. Introduction A submerged ¯oating tunnel can be an alternative crossing both to rock tunnels below the sea bottom and to long span bridges above the surface. Compared with long span cable-stayed and suspension bridges, the increase of cost with span length is generally much better for the submerged ¯oating tunnels. Rock tunnels will on the other hand increase in cost for necessary depth to rock bottom (increasing the total length) and the rock conditions with respect to excavation costs. The submerged tunnel will also be favorable if interference with the landscape and scenic view is to be minimized. The submerged ¯oating tunnel in Fig. 1 is proposed for the crossing of Hùgsfjorden in Rogaland County, Norway. The position of the tunnel relative to the surface is for this concept maintained by buoyancy forces provided by the pontoons. Alternatively the tunnel can

* Corresponding author.

be moored to the bottom, or kept in position by a combination of mooring lines and pontoons. The alternative thought to be the most feasible, will depend on the site conditions such as span length, water depth, bottom conditions, ship trac, environmental loading (waves, current, earthquake) etc. The ®rst feasibility study for a submerged ¯oating tunnel in Norway was reported to the Norwegian Public Roads Administration in 1971 for a speci®c fjord crossing. The concept was based on cable mooring to the bottom at two locations along the span. The submerged tunnels in Figs. 1 and 12 were proposed in response to the pilot study for Hùgsfjorden in 1988± 1989. The present paper discusses the merits of various methods for stochastic dynamic response analysis of a submerged ¯oating tunnel subjected to wind driven surface waves. The choice of method will depend on the actual structural concept but also on the particular features required for the di€erent stages in the design process. Di€erent levels of approximations will be rational for a concept study vs the ®nal design analyses

0045-7949/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 3 2 9 - 0

660

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 1. Possible concept for a submerged tunnel across Hùgsfjorden.

to check if the speci®ed safety criteria are met. The highest level of method development is represented by the Navier±Stokes formulation as implemented in a ®nite element (FE) analysis of the ¯uid/structure interaction (FSI). Such computational ¯uid dynamics (CFD) procedures are expected to provide reference solutions for ¯uid/structure interactions of submerged ¯oating tunnels in the near future. Particular interest is associated with components with typical three-dimensional (3-D) ¯ow. Examples are pontoons and connections between the tunnel and the pontoons, see Fig. 1. Such reference solutions will contribute signi®cantly to obtain more economic structures still meeting the target safety levels. When Navier±Stokes FSI analyses become robust and ecient, a complete package of computer analysis tools with a sound theoretical basis will be available. This will complete a continuous development and improvement of methods of analysis over a couple of decades. The numerical examples in the present study focus on some global dynamic response characteristics of submerged ¯oating tunnels exposed to wind driven waves. Emphasis is also put on the performance of frequency domain vs time domain formulations for the stochastic dynamic analysis. The choice between the alternatives will depend on the type of structural concept. In addition, computation with the Navier±Stokes FE model is demonstrated for a two-dimensional (2-D)

sectional model of the tunnel cross section in regular waves. 2. Stochastic dynamic response analysis 2.1. Dynamic equilibrium equation The relationship between the load and response stochastic processes may be expressed in terms of the semi-discrete equations of motion: RI ‡ RD ‡ RS ˆ Q I

…2:1† D

R is the inertia force vector, R is the damping force vector, RS is the internal force vector and Q is external force vector. The above relation states the instantaneous equation of motion of a non-linear system at time t. The inertia forces may be expressed in terms of the mass matrix, M, and the accelerations, rÈ : … ‡1 RI ˆ RI …Èr ,t† ˆ M…t ÿ t†Èr dt …2:2† 0

The work done by the viscous damping force vector, RD, represents the total energy dissipation in the system. It is established by considering the tangential damping matrix CT. The equivalent damping force may be written as:

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

RD ˆ RD …Èr ,Çr ,r,t† ˆ

… ‡1 0

CT …Çr ,r,t ÿ t†Çr dt

…2:3†

The lower integration limit in the expressions above is zero, due to the principle of causality. The convolution integrals are due to the mass and damping matrices being frequency dependent. Non-linearity is accounted for by allowing the damping matrix to be a function of displacements and velocities, the internal force vector, a function of displacements and the external force vector, a function of displacements, velocities and accelerations, respectively. The mass matrix includes both structural and hydrodynamic contributions, with the latter being frequency dependent while the former is generally taken to be constant for submerged tunnels. Similarly, the damping matrix includes a frequency dependent hydrodynamic term and a constant structural contribution modeled from estimates of modal structural modal damping ratios. Typically, the dependence of the displacement and velocity vector is neglected for the present type of structure. The structural sti€ness matrix also includes contributions due to possible surface-piercing members such as pontoons in addition to tethers and anchor lines. For the special case that the frequency dependence of the mass and damping matrices can be neglected, the two convolution integrals in the equilibrium equation disappear, resulting in a simple product form for the inertia and damping terms. For a stepwise time-integration of the equations of motion, it is convenient to express them on their corresponding incremental form: MDÈr …t† ‡ CT DÇr …t† ‡ KT Dr…t† ˆ DQ…t†

…2:4†

where M, CT and KT are the mass, the tangential damping and the tangential sti€ness matrices valid for each incremental step; the D symbol designates increment of accelerations, velocities, displacements and load vectors over the time step. This equation will also represent non-linear response due to geometric sti€ness. Techniques for ecient solution of the equations of motion are frequently based on schemes for a reduction of the dimension of the equation system. One such method is transformation into normal coordinates. However, in the numerical examples to follow, the full size of the equation system is kept. 2.2. Response analysis in the frequency domain 2.2.1. Relationship between load and response spectral matrices In the following, it is assumed that the structural matrices and the hydrodynamic load has been linearized.

661

From the theory of multi-degree-of-freedom systems under random loading, the response spectral density matrix can be expressed in terms of the load spectral matrix as: Sr …o † ˆ H…o †SQ HT …o †

…2:5†

here H…o † ˆ ‰KT ‡ io CT …o † ÿ o 2 M…o †Šÿ1

…2:6†

is the complex frequency response function of the structure-¯uid system. In turn, the load matrix can be expressed by the one-dimensional sea elevation spectral density, yielding: Sr …o † ˆ H…o †F…o †HT …o †SZ …o † ˆ B…o †SZ …o †

…2:7†

where F(o ) is referred to as the hydrodynamic transfer function. The matrix B(o ) is designated the transfer function of the total system relating the response spectral density matrix to the scalar wave spectral density. These two matrices also account for directionality e€ects related to hydrodynamic loading and structural response. In addition to the mean direction of the incoming wave system, the directional spreading of wave energy around this direction is important. This is expressed in terms of the wave spreading function C(y ) where y is the wave direction. The load amplitudes (per unit wave height) will, typically, also depend on the direction of the incoming wave. Furthermore, the structural response may exhibit a strong sensitivity to mean wave load direction. By introducing the direction and frequency-dependent loading Q(o, y ), the directional integration can be expressed on the following form: … F…o † ˆ Q…o ,y†QT …o ,y†C…y† dy …2:8† y

where the directional integration can be performed either on element or system level. Accurate computation of the direction and frequency dependent load vector in general requires utilization of advanced and time-consuming numerical schemes such as the sink-and-source method. The same applies to the frequency-dependent hydrodynamic contributions to the damping and mass matrices. In particular, submerged tunnels with complex geometries, caused, for example, by surface pontoons, demand careful modeling and ®ne-meshed spatial discretization if high precision is aimed at. However, at a preliminary design stage simpli®cations can be introduced in order to increase the computational eciency. As an example, three-dimensional hydrodynamic models can be assembled from a sequence of twodimensional approximations. Furthermore, analytical

662

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

or semi-analytical expressions for axi-symmetric threedimensional components such as pontoons and shafts can be utilized. 2.2.2. Response statistics for design purposes Due to the linear relationship between the load and response stochastic processes, the Gaussian statistical properties assumed for the former is also preserved for the latter process. Extreme value statistics for Gaussian processes are typically modeled by a Gumbel distribution. The statistical parameters of this distribution are functions of the expected number of local response maxima which occur during a given reference time period (e.g. an extreme wave condition with a duration of some hours). This number is expressed as: Nˆ

s_ 2ps

…2:9†

where s is the standard deviation of the response process, and the superscript dot in the numerator refers to response velocity. The expected value of the Gumbel distribution (i.e. the expected largest response for the given duration) is expressed as:  p 0:5772 E‰x max Š ˆ s 2 ln N ‡ p 2 ln N

…2:10†

However, in some cases the most probable largest response is employed, which corresponds to keeping only the ®rst term in the parentheses of this expression. The spread around the expected value is expressed by the standard deviation: p s‰x max Š ˆ s p 2 3 ln N

…2:11†

which decreases for increasing values of N. Variances and corresponding standard deviations of response displacements are obtained directly from the spectral response matrix by integrating the proper response spectral density with respect to frequency. The spectral densities for the cross-section forces and stresses are derived from these by pre- and post-multiplication with the proper element sti€ness matrices and cross-sectional resistance moments. Due to the dynamic ampli®cation e€ects, there will be a signi®cant sensitivity of the response with respect to a variation in a peak period of the extreme sea state. Accordingly, a more consistent approach is to establish a so-called long-term statistical distribution. Probability distributions of local response maxima is then computed for a number of sea-states. The signi®cant wave heights and peak periods are identi®ed, together with corresponding frequencies of occurrence, from the so-called scatter diagram applicable for the

relevant bridge site. The long-term distribution of local response maxima is subsequently obtained as a weighted combination of the distributions for all these sea-states. The response level corresponding to a given return period is ®nally estimated by specifying the corresponding probability of exceedance and inverting the long-term distribution. A particular feature of some submerged ¯oating tunnel bridge concepts is also a signi®cant sensitivity of structural response to mean wave direction and degree of directional spreading, as illustrated below. This implies that the probability of occurrence for a number of wave directions and corresponding spreading parameters must be considered when computing the longterm distribution. 2.3. Response analysis in the time domain 2.3.1. Computation of response time series The time domain analysis is based on simulation of one or more samples of the stochastic load process. The samples can be generated, for example, by means of Monte Carlo simulation methods, see Refs. [1±3]. For each of these sample functions, the corresponding load vector time series is computed, and subsequently a deterministic type of response analysis is performed. In this type of approach, the sea-elevation process can be approximated by a discrete sum as: Z…x,t† ˆ

N1 X N2 X

Aml fcos o m t

mˆ1 lˆ1

…2:12†

ÿ k…o m †…x cos yl ‡ y sin yl † ‡ fml g where Aml ˆ

p 2SZ …o m ,yl †Do Dy

…2:13†

Here, jml are random independent phase angles distributed uniformly between 0 and 2p, and k is the wave-number. Furthermore, Do ˆ

Dy ˆ

o up N1

and o m ˆ …m ÿ l †Do

yup ÿ ylow N2

and

yl ˆ ylow ‡ …l ÿ 1†Dy

…2:14a†

…2:14b†

where oup is the upper limit of the frequencies to be included in the analysis; yup and ylow designate the upper and lower limits for the summation with respect to wave direction. The double summation can be eciently carried out by means of FFT-techniques. Evaluation of the load vector for the submerged bridge can be performed

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

based on the same type of expression simply by inserting the frequency and direction dependent (complexvalued) vector transfer functions between each harmonic component of the sea-elevation process and the corresponding harmonic load vector component, i.e.: Q…t† ˆ

N1 X N2 X

Aml Q…o m ,yl †fcos o m t

mˆ1 lˆ1

where for each component of the load vector, a double summation of the same type as for the sea-elevation process is obtained. For each generated sample of hydrodynamic load time series, corresponding response time series are computed by incremental step-by-step integration. A general incremental-iterative scheme combined with the Hilber±Hughes±Taylor time-approximate recursion algorithm [4] is used for solving the equations of motion discretized in space. Vectors of incremental accelerations and velocities are expressed by the vector of incremental displacements and the accelerations and the velocities at the previous time step, such that: Dr n‡1 ˆ r n‡1 ÿ r n ˆ

1  1 r_ n ÿ r n Dt2 b 2b

Dr_ n‡1 ˆ r_ n‡1 ÿ r_ n ˆ

g g Dr n‡1 ÿ r_ n ÿ D Dtb b

t



 g ÿ 1 r n 2b

…2:16†

…2:17†

n is the time step number, Dt=tn+1ÿtn the time increment, and g and b are parameters de®ned by the Hilber±Hughes±Taylor parameter, aH, as 1 g ˆ …1 ÿ 2aH † 2

…2:18a†

1 …1 ÿ aH †2 4

…2:18b†

bˆ and

1 ÿ
DQn‡1 ˆ …1 ‡ aH †‰Qn‡1 ÿ Qn ‡ CT,n bn Š ‡ Mn an Än ‡Q

…2:18c†

rn, Çrn and rÈ n are the approximations to r(tn ), Çr(tn ) and rÈ (tn ), respectively. The associated incremental load vector is found to be

…2:19†

where an ˆ

1  r_ n ‡ Dtb

 1 ÿ 1 r n 2b

…2:20a†

bn ˆ

  g g r_ n ‡ Dt ÿ 1 r n b 2b

…2:20b†

…2:15†

ÿ k…o m †…x cos yl ‡ y sin yl † ‡ fml g

663



Ä n ˆ Qn ÿ RS ÿ CT,n r_ n Q n

…2:20c†

Qn and RSn are external and internal force vectors, while Mn and CT,n are the mass and the tangential damping matrices at step n. Eq. (2.19) accounts for the force unbalance at time step n, in addition to the increment of external loads for time step (n+1). Unbalance in the equations of motion (dynamic equilibrium) will thus not be accumulated. 2.3.2. Estimation of response statistics for design purposes For each sample response time series (obtained for a sea-state with given signi®cant wave height, peak period, wave direction and wave spreading), one single value representing the largest response for the given duration is identi®ed. However, based on all the given local response maxima of the time series, corresponding probability density and distribution functions can be estimated, for example, by the ®tting of proper theoretical models to the empirical data. The parameters of these ®tted models can subsequently be employed for identi®cation of relevant extreme value distributions. Furthermore, long-term distributions of relevant response quantities are obtained by weighting the local maxima distributions for a number of seastates in a manner identical to the frequency domain procedure. However, for the time domain analysis, the required simulation length to achieve adequate con®dence in the estimated statistical parameters becomes a critical issue. This is of particular importance for cases where non-linear behavior is expected to be important. For such cases, estimation of the upper fractiles of the local maxima and extreme value distributions is very sensitive to the number of response samples available. Generally, this implies that a number of sample response time series are generated for each sea-state. The mean values of the statistical parameters obtained from the di€erent samples are then employed. At the same time, the spread of these parameters between

664

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

samples allows an assessment of the required number of repetitions to achieve a desired accuracy. 2.4. Assessment of response analysis procedures in the frequency vs time domain

submerged tunnel. The objective of these simulations is to obtain inertia, damping and excitation force coecients for the whole cross section, which are utilized in global analyses for structural parts exposed to 3-D ¯ow.

Summarizing the main features of the response analysis procedures, the bene®ts of the frequency domain approach are as follows: . frequency dependency of hydrodynamic coecients for damping, added mass and excitation forces are easily incorporated; . estimation of extreme values for design purposes is straightforward due to available closed-form analytical expressions; . the computational e€ort is signi®cantly reduced as compared to time domain analysis. However, on the negative side, this type of analysis is not well suited for incorporation of non-linearities related to structural behavior or hydrodynamic modeling. The bene®ts of the time domain approach can be summarized as: . non-linear e€ects related, for example, to non-linear material behavior, geometric sti€ness, ®nite surface wave e€ects and viscous loading can be incorporated in a direct manner; . simulation in the time domain provides insight into the physical behavior of the structure.Instantaneous deformation patterns and time variation of response can be easily visualized. Among the negative aspects of this procedure are; increased e€ort related to estimation of extreme response statistics, and diculties related to implementation of frequency dependent damping and mass coecients. Typically, both types of analyses should be performed in order to assess the importance of various modeling assumptions in relation to a computed response. 3. Navier±Stokes modeling of ¯uid/structure interaction True behaviour of the ¯oating tunnel problem can be captured by utilizing a computational tool that combines the ¯uid, the structure and their interaction in an integrated single-pass simulation. In the current study, a 2-D sectional model of the ¯oating tunnel problem is analyzed with the multiphysics FE code SPECTRUM [5]. The model consists of a circular cylinder submerged in sea water and subjected to constant current and regular waves. Structural boundary conditions are consistent with a global model of the

3.1. Computational model The case considered herein corresponds to a straight tubular bridge with a diameter of 10 m and a length of 800 m. The tube is supported at the two ends only (no pontoons or tension legs). The conditions simulate a water depth of 100 m with the center of the tube being 30 m below the water surface. The FE model of this coupled FSI-problem consists basically of three main components (1) a ¯uid region where the incompressible Navier±Stokes equations are to be solved; (2) a solid region where (geometrically) non-linear elasticity equations are to be solved; (3) an interface de®ning the coupling between the two regions. A complicating factor in the solution process is that the boundary between the two regions is not ®xed in space; it moves according to the total displacement of the structure. In addition, the ¯uid boundary representing the water surface must be permitted to move in order to represent the wave motion properly. This necessitates the Arbitrary Lagrangian Euler (ALE) formulation for the ¯uid analysis [6], where the FE mesh is permitted to move independently of the motion of the ¯uid particles themselves.

3.1.1. The ¯uid region The geometry of the ¯uid region is shown in Fig. 2. The problem is solved here as a 2-D problem using one hexahedron element in the length direction of the cylinder (global y-direction). The thickness of the ¯uid domain is set to L=10.0 m. The ¯uid is assumed to be sea water with density rw=1025.0 kg/m3 and dynamic viscosity m=0.0016 kg/ms. A constant vertical body force is applied through the gravity constant g=9.81 m/s2. The mesh used in the current study, is shown in Fig. 3. The mesh is graded towards the cylinder where a boundary layer develops due to viscous e€ects. 3.1.1.1. The incompressible Navier±Stokes equations. The ®eld equations for which the solution is sought in the ¯uid domain are as follows:

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

665

Fig. 2. Geometry of the ¯uid region. All units in metres.

rw uÇ ‡ rw …u
r†u ÿ r
s ÿ r
p ˆ f

…3:1†

r
u ˆ 0

…3:2†

where u={ux, uz }T (the velocity vector) and p (the pressure) are the unknown variables, f={0, ÿrwg }T is the speci®ed body-force vector and   @ @ ix , iz rˆ @x @z denotes the gradient operator. The components of the viscous stress tensor, s, are given as sij ˆ m…ui,j ‡ uj,i †,

i,j 2 fx,zg

…3:3†

In addition, the ALE mesh displacements (d ALE and x d ALE ) are solved from a separate equation of large z

deformation elasticity. There is therefore a total of ®ve unknown variables (for the 2-D sectional model) at each node of the FE mesh. The equations above are solved in space and time in SPECTRUM, using an implicit FE method and Newmark time integration. A system of non-linear algebraic equations then arise for each time step, which are solved by the Broyden, Fletcher, Goldfarb and Shanno (BFGS) quasi-Newton method. 3.1.1.2. Boundary and initial conditions. The following nodal boundary conditions (Dirichlet or essential conditions) are applied on the four ¯uid boundaries: 1. in¯ow boundary ux=uc+uw(z, t ), uz=ww(z, t ) and d ALE =0; x 2. out¯ow boundary and free surface boundary:

Fig. 3. Discretizations of the ¯uid region: FE mesh with 832 elements and 2  893 nodes.

666

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

d ALE =0; x 3. bottom boundary: uz=0, d ALE =0 and d ALE =0; x z 4. wave excitation nodes: (i.e. at x=z=0): d ALE =d w z z (z, t ). Here, uc is the (constant) current velocity, uw(z, t ) and ww(z, t ) denote the ¯uid particle velocity at the in¯ow boundary according to linear wave theory and d w z (z, t ) is the corresponding vertical displacement of the ¯uid particles. These three functions are given by Ref. [7] uw …z,t† ˆ o za

cosh k…z ‡ h† sin o t sinh kh

…3:4a†

sinh k…z ‡ h† cos o t sinh kh

…3:4b†

ww …z,t† ˆ o za

dw z …z,t† ˆ za

sinh k…z ‡ h† sin o t sinh kh

…3:4c†

where za is the wave amplitude, o=2p/T is the angular frequency (T is the wave period), h is the water depth and k is the wave-number. The latter is determined from the equation o 2 ˆ gk tanh kh

…3:5†

It is emphasized that the horizontal mesh movement, d ALE , is prescribed equal to zero along all outer x boundaries. This is, however, not inconsistent with the non-zero horizontal ¯uid velocities, ux, since in the ALE formulation there is no direct coupling between these two quantities. The following element boundary conditions (Neumann or natural conditions) are applied:

1. out¯ow boundary: pressure p and tangential tractions are `internally de®ned'; =un. 2. free surface boundary: p=0 and u ALE n The notion `internally de®ned' means that the value of the boundary condition in question is to be computed from the solution at the boundary and may be regarded as a `no boundary condition' condition. The free surface condition u ALE =un states that the normal n velocity of the mesh is set equal to the normal component of the ¯uid velocity, i.e. no ¯uid is allowed to ¯ow through the moving boundary. The free surface (where p 0 0) and the gravity load give a linearly varying pressure equal to p=ÿrwgz (hydrostatic pressure) when the velocity is zero. This will also be the dominating part of the total pressure when the velocity is non-zero. The hydrostatic pressure is therefore used as initial condition on the pressure in the simulations. For the velocity, we use ux=uc and uz=0 as initial conditions. 3.1.2. The solid region The structural part of the FSI-problem consists of a straight beam in the y-direction which is completely ®xed at the two ends. The beam has a circular hollow cross section with outer diameter D=10.0 m and wall thickness t=1.0 m, and is neutrally buoyant in water. A hypoelastic stress model is assumed for the cylinder wall with material parameters E=26.9  109 N/m2 and n=0.2. In a 2-D simulation of the coupled FSI-problem the beam structure is ideally modeled as a simple two degrees-of-freedom (DOFs) system with a point mass M=rwAwL (where Aw and L are the cross sectional area and length of the displaced ¯uid body, respect-

Fig. 4. Geometry and properties of the solid region. (a) Idealized model. (b) The model used in the present SPECTRUM simulations.

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

ively), and where the global bending sti€ness of the beam is represented by two identical springs, see Fig. 4(a). However, in order to de®ne the interface between the ¯uid and solid region in SPECTRUM, a section of the beam must be modeled using either solid or shell elements. The program does not o€er the possibility to de®ne an interface between ¯uid continuum elements and structural beam elements, so the beam section cannot be modeled using a single beam element (which perhaps would have been preferable for this case). In the present study, the beam cross section is modeled using 16 solid elements as shown in Fig. 4(b). These elements are then attached to a system of truss elements that represents the global bending sti€ness of the beam. To prevent the circular beam section from spinning around its own axis, the vertical spring is divided into two truss elements attached to opposite sides of the cylinder. Note that since we neither use shell elements nor beam elements we avoid rotational DOFs in the FE model. The sti€ness and mass of the solid model, represented by the spring sti€ness ks and the solid density rs, respectively, must be selected such that the overall static and dynamic properties of the 2-D sectional model correspond to those of the global structure. We ®rst determine the spring sti€ness from a static consideration. This sti€ness is simply the ratio between the total load qL acting on the beam segment in the direction of the spring (q denotes the load intensity and L is the length of the segment) and the de¯ection d of the beam section. For a segment located at the midpoint of a beam of total length l a pure linear analysis then yields the sti€ness

ks ˆ

qL EI ˆ 384 4 L d l

…3:6†

where I=p/64[D 4ÿ(Dÿ2t )4] is the moment of inertia of the beam cross section. For a segment of length L=10.0 m, of a beam with total length l=800.0 m, this yields ks=73,087 N/m, with the cross sectional data given above. However, since the beam is ®xed at both ends also in the length direction, a transverse de¯ection will introduce a tension force in the beam resulting in an additional geometrical sti€ness. Therefore, to obtain a more accurate estimate of the true sti€ness, a geometrically non-linear analysis of the global structure using the FE code FENRIS [8] is performed. This non-linear analysis results in ks=82,552 N/m for de¯ections of the order of d 1 1.0 m. The total e€ective mass M of the solid region is now determined such that its natural frequency corresponds to the lowest frequency, o1 of the global structure, i.e.:

ks o 21

Mˆ

667

…3:7†

We here assume an added mass coecient Cm=1.0 such that the e€ective mass is given by M ˆ ‰rs As ‡ …1 ‡ Cm †rw Aw ŠL

…3:8†

where As denotes the cross sectional area of the solid domain Aw and is the cross sectional area of the de¯ected ¯uid body. By combining the two Eqs. (3.7) and (3.8) we now obtain the following expression for the solid density rs ˆ

ks =o 21 ÿ …1 ‡ Cm †rw Aw L As L

…3:9†

Note that since the circular cross section here is discretized with linear elements, the actual areas As and Aw are slightly smaller than the area of a perfect circle. If ncs and ncw denote the number of (equally sized) elements around the circle for the solid and ¯uid domains, respectively, these two areas are given by, respectively As ˆ …D ÿ t†tncs cos

Aw ˆ

p p sin nns ncs

D2 p p ncw cos sin ncw ncw 4

…3:10†

…3:11†

An eigenvalue analysis of the submerged global structure is now performed using FENRIS, resulting in the lowest angular frequency o1=0.138 rad/s. With ncs=16, ncw=32 and L=10.0 m, Eqs. (3.9), (3.10) and (3.11) then yield the solid density rs=9884.5 kg/m3. Finally, we determine a vertical body force, bz, for the solid region such that the cylinder is neutrally buoyant, i.e.: bz ˆ g

rw Aw ˆ ÿ2:88 m=s2 rs As

…3:12†

3.1.3. Fluid-structure interface The DOFs of the two regions are tied together by means of an interface de®ned on the cylindrical boundary between the two regions. The interface is de®ned by specifying the nodal connectivities of the two boundary surfaces of the respective regions. The boundary surface of the solid region is de®ned as the master surface, i.e. the surface for which contact/penetration of nodes belonging to the other (slave) surface is checked. The interface conditions are then formulated as follows: 1. the nodes on the slave (¯uid) surface are tied in

668

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

®xed positions with respect to the master (solid) surface regardless of voids between the two regions due to di€erent levels of discretization; 2. no-slip velocity conditions, i.e. u¯uid=usolid on the interface; 3. no-slip mesh conditions, i.e. dALE=dsolid on the interface. The velocity and mesh conditions may be enforced either through Lagrangian multipliers or by penalty factors. In the present study the penalty factor approach is used with the penalty factors 100.0 and 1.0 for the velocity and mesh constraints, respectively. 3.2. Solution strategy The solution strategy in the SPECTRUM solver consists of several nested loops, with a number of control parameters to be speci®ed by the analyst at each level, see Fig. 5. Although the program provides default choices for each parameter it is necessary to ®ne-tune several parameters in order to obtain proper results. The amount of tuning depends on the complexity of the problem that is analyzed. The particular problem considered in the current work (involving both FSI and free surface phenomena) turned out to represent a challenge for this code.

After thorough experimentation, a solution strategy that works fairly well for the 2-D sectional model of the ¯oating tunnel problem was established. The key features of this strategy are as follows: . Two staggers are employed (i.e. two equation subsystems where the corresponding ®eld equation is solved) Stagger 1: Incompressible isothermal Navier±Stokes equation in the ¯uid region (i.e. primary variables ux, uz, and p ) and linear isotropic elasticity equation (in®nitesimal strains) in the solid region (i.e. primary variables dx and dz ). Stagger 2: ALE mesh movement equation in the ¯uid region (i.e. primary variables d ALE and d ALE ). x z When a particular stagger is solved, the active variables of the other stagger are `frozen' and act as input parameters to the active stagger. Once converged, the solution variables are updated and are then passed to the other stagger as frozen variables. This process is then repeated sequentially within the Stagger Control Loop until both staggers have converged; . The convergence tolerance is speci®ed independently for each equation of the problem, in terms of the residual norm and/or in terms of the solution increment norm. In the current work, the following tolerances are used: Equation Navier±Stokes Elasticity Mesh movement

Fig. 5. Schematic overview of the SPECTRUM solution procedure.

Residual 0.001 0.001 none

Solution increment 0.01 0.001 0.001

. Stagger 1: One non-linear iteration and maximum two augmented Lagrangian iterations. The lefthand-side matrix is recalculated for each iteration. Note that in the non-linear loop (the innermost loop in Fig. 5) only the unconstrained DOFs are considered. The nodal boundary conditions (and the interface conditions unless the penalty factor approach is used for those) are enforced using augmented Lagrangian multipliers which are updated only in the outer augmented Lagrangian loop. . Stagger 2: Maximum ten non-linear iterations with an update of the left-hand-side matrix only in the ®rst iteration, and one augmented Lagrangian iteration. This stagger converges usually very fast (within 2±3 iterations). . Stagger Control Loop: Maximum ten iterations. However, most time increments converge within 4-5 iterations in the present simulations. . Linear solver: A sparse direct solver is used for both

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

669

staggers. An iterative solver that may be run in parallel by means of domain decomposition is also available. However, the size of the FE model must be larger than that used herein for that solver to be competitive. The main characteristics of the above strategy is that very few iterations are performed internally for one stagger before moving to the next stagger. Instead we allow for more stagger iterations (i.e. iterations in the outer Stagger Control Loop). Consequently, each stagger usually have not fully converged when we switch to the other stagger, except for in the ®nal stagger iteration. This strategy results in a limited overall number of iterations for each time increment and is therefore found to be the more cost ecient. 4. Numerical studies 4.1. Global response from wind generated waves Loading due to wind generated waves is a major load case for the assessment of dynamic response properties of a submerged ¯oating tunnel. The following numerical studies will focus on the e€ect of spatial variation of the loading (mean wave direction and wave spreading) and in¯uence on estimated response from hydrodynamic mass, damping and geometric sti€ness. 4.1.1. Directionality e€ects for straight bridge and horizontally curved bridge The straight bridge is chosen for illustration of the principal e€ect. Similar e€ects are then illustrated for curved bridges which represent the most likely geometries for a real submerged tunnel. The straight bridge structure is modeled by 20 3-D beam elements of equal length, which is satisfactory with respect to representation of the natural modes of vibration contributing to the wave response. Also the load correlation along the bridge is well represented within this structural model. The response analyses are performed in the frequency domain with the computer program FEDAF [9]. The hydrodynamic loading and ¯uid interaction is calculated on the basis of linear potential theory. The JONSWAP [10] model spectrum has been used for the sea elevation. The directional spreading of wave energy is modeled by a frequency independent cosine function. Two of the natural modes of vibration in the vertical direction for the straight bridge are illustrated in Fig. 6. The corresponding natural frequencies are o7=1.11 rad/s and o8=1.287 rad/s. We apply long crested wave loading represented by an elevation spectrum with peak circular frequency op=0.4p (1.26 rad/s). Based

Fig. 6. Natural modes for straight bridge.

on the frequency ratios, signi®cant dynamic ampli®cation is anticipated for both modes. The degree of dynamic ampli®cation is, however, also related to the spatial load variation relative to the mode shapes. Maximum `spatial resonance' is obtained when the wavelength of a given mode coincides with the wavelength of the incoming wave decomposed onto the bridge axis. This is illustrated in Fig. 7 for

Fig. 7. Parallel wave crests propagating at an angle y from the bridge normal.

670

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 8. Maximum expectation value of bending moment in vertical bridge motion as function of mean wave direction.

long crested waves propagating in a direction de®ned by the angle y. The cross sections of the bridge at a distance Ly apart will be in phase with respect to hydrodynamic force while the mid point of the Ly portion of the bridge will be 1808 out of phase. In order to obtain maximum dynamic ampli®cation in mode shape 8 the value of Ly should be equal to the distance between the outer peaks of the mode. From Fig. 6 this is esti-

Fig. 10. Modal load vector amplitudes versus wave propagation direction.

mated as Ly=480 m, corresponding to y=4.678 (see Fig. 7) for a wavelength lp=39 m (lp=0.4p). For mode 7 the corresponding Ly will be twice the distance between the peaks, giving Ly=700 m and y=3.28. Around y=48 we should then expect to get considerable dynamic ampli®cation from both modes 7 and 8. This is veri®ed by numerical calculations as showed in

Fig. 9. Expected peak response along the bridge for di€erent wave propagation directions.

S. Remseth et al. / Computers and Structures 72 (1999) 659±685 Table 1 Possible wave propagation angles to give ampli®ed response in modes 7 and 8 ModeLy (m)y1 (8) 12 Ly (m)y1/2 (8) 14 Ly (m)y1/4 (8) 16 Ly (m)y1/6 (8) 7 8

700 480

3.2 350 4.67 240

6.4 9.35

175 120

12.9 19.0

116.7 80.0

19.5 29.2

Fig. 8, where the maximum response along the bridge axis has been selected for a given y. If we look at the distribution of expected peak response along the bridge we can clearly identify from Fig. 9 the change from the symmetric and low level response of 08 towards a much higher level and asymmetric response at 38, 48 and 58 angle of wave propagation. The shape of the curves of expected peak response for the latter angles is also easily identi®ed as a combination of modes 7 and 8. A qualitative veri®cation of the numerical results can be obtained analytically by considering an idealized case. A modal analysis of a simply supported beam with mode shapes consisting of an increasing number of half sine waves is easily performed by the Rayleigh±Ritz technique. Results from these calculations are shown in Fig. 10 con®rming the said directional sensitivity for a straight bridge. Some ampli®cation of modes 7 and 8 can be expected also for angles y corresponding to 12 Ly, 14 Ly, 1 6 Ly and so on, especially if the possibility of an interaction between the two modes exists such as for y 1 48. Comparing Table 1 with results shown in Fig. 8, we may identify an ampli®cation close to y 1 108 corresponding to 12 Ly for mode 8. A plot of expected peak response along the bridge is also symmetric for this angle. Furthermore, we observe a distinct peak in Fig. 8 for y 1 208. The response distribution along the bridge, then has approximately the same form as for y

Fig. 11. Maximum expectation value of bending moment about horizontal axis versus wave direction.

671

1 48 (see Fig. 9). From Table 1 we obtain 14 Ly for mode 8 and 16 Ly for mode 7 at approximately 208. Interaction between the modes can thus be expected and explains the second peak in Fig. 8. For a curved bridge with the same distance between the end points but with a radius of curvature equal to 1700 m, the maximum bending moments about a horizontal cross sectional axis as function of wave propagation angle is shown in Fig. 11. It is seen from this ®gure that we get signi®cant e€ects from mean wave direction in the vicinity of the normal direction for a curved bridge. The maximum bending moment is, however, lower for this case than for the straight bridge. Due to phase e€ects, the modal loads are expected to decrease. For the curved bridge there is also a combined vertical/horizontal symmetric mode shape with a frequency ratio equal to 0.9. This seems to contribute to give the maximum expected bending moment midway between the end points for y=0, and the response for y=0 is also three times as large as for the straight bridge for this angle. The second peak in Fig. 11 indicates that also curved bridges may require investigation of mean wave direction in the vicinity of the normal direction. If we look at the bending moment about a vertical axis this is con®rmed as the maximum there will occur for y 1 68. 4.1.2. E€ects of wave direction and shortcrestedness for a vertically curved tunnel with tether mooring The vertically curved bridge moored with tension legs to anchors on the bottom is one of the possible concepts for crossing of Hùgsfjorden, see Figs. 1 and 12. The overall spanwidth is 1500 m and the bridge has four equally spaced tension legs. The net buoyancy of the tunnel is approximately 10% of the submerged weight. The bridge is modeled by 50 equally sized 3-D beam elements. The geometric sti€ness from the tension legs are modeled by linear elastic springs with constant sti€ness. The analyses are performed in the frequency domain with the computer program FEDAF. Figs. 13 and 14 show a couple of natural mode shapes for the bridge in dominant sway and heave motion respectively. The wind driven wave loading is modeled with a JONSWAP sea elevation spectrum with peak period Tp=5.5 s. For long crested waves at a zero angle of incidence, considerable dynamic ampli®cation of the wave response should be expected for mode 5 in heave and mode 6 in sway. These modes are symmetric and with initial frequency ratios at Tp of 1.04 and 0.96 respectively. The damping will also be light in these modes with the structural ratio set to 0.4% and limited displacement response (less than 0.05 diameter) for hydrodynamic energy dissipation. The in¯uence of the angle of incidence on the bend-

672

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 12. Submerged ¯oating tunnel curved in a verticle plane and with tension leg mooring.

ing moment response is discussed for sway motion where mode 6 is supposed to contribute considerably to the total response according to the discussion above. For a zero angle of incidence the expected peak bending moment response along the bridge due to sway motion is shown in Fig. 15. It is clearly seen by comparison with mode shape 6 that apart from the ®xed end moments, the maximum values correspond to the shape of this mode. (Note that the de®nition of expected bending moment response is always positive.) The bending moment magnitude of the largest peak (excluding the ®xed ends) is approximately 1.85  106

kNm. Fig. 16 gives equivalent results for an angle of incidence of 18 the normal. The largest peak is now at 0.5  106 kNm. However, it is noticed that the locations of the peaks correspond with the asymmetrical mode shape no. 4. In order to explain the large decrease in response with the small change of angle of incidence, there are two major factors to consider. One is the Ly length along the bridge (see Section 4.1.1) between equal wave loading conditions from two consecutive waves of non-normal direction to the plane of the bridge. The other factor is the variation along the tunnel of the distance between the surface and the top of the tunnel cross section. At the ends this distance is

Fig. 13. Natural modes of vibration in heave for the tunel in Fig. 12.

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

673

Fig. 14. Natural modes of vibration in sway for the tunnel in Fig. 12

7 m and with a variation de®ned by a circular arch down to 25 m at the mid section, see Fig. 12. The wave length corresponding to the peak period will be 31.5 m. This will hardly give any wave loading in the

middle of the bridge, but with increasing magnitude toward the ends. The value of Ly for 18 from the normal direction will be 1805 m for l=31.5 m. This will provide no extra dynamic ampli®cation of the sym-

Fig. 15. Expected peak level of bending moments from sway motion due to long crested waves at perpendicular incidence.

674

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 16. Expected peak level of bending moments from sway motion due to long crested waves at 18 from perpendicular incidence.

metric mode shapes. At 1/2Ly on the other hand, the loading from two consecutive waves will be 1808 out of phase and can thus amplify the response in an asymmetrical mode shape. The distance of 1/2Ly 1 900 m may be quite close to the distance between load

resultants in opposite directions at the two ends of the tunnel. The corresponding distance between the outer peaks of mode 4 is 1000 m. Excitation of the asymmetric mode is thus likely. The frequency ratio at Tp is 1.4 for mode 4. The total response is a combination

Fig. 17. Expected peak level of bending moments from sway motion due to long crested waves at 1.58 from perpendicular incidence.

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

675

Table 2 Added mass, wave damping, structural damping, modeling of tension legs and net buoyancy for the di€erent dynamic analysis carried out (T; coecients picked from Tables in Ref. [11] for a frequency corresponding to the peak period Tp=5.5 s, C; tension legs modeled by non-linear cable elements, S; tension legs modeled by linear springs) Case

Cm

B

aM=aK

Tension legs

Buoyancy (%)

1 2 3 4 5 6 7 8

T T 1.0 1.0 T T T T

0.0 T 0.0 T 0.0 T T T

0.004 0.004 0.004 0.004 0.01 0.01 0.004 0.004

C C C C C C S C

10 10 10 10 10 10 10 5

with symmetrical mode shapes, but according to the expected peak response in Fig. 16 the asymmetric load e€ect is dominating for the bending moment from sway motion. Fig. 17 gives the expected peak response of the bending moment due to sway motion for an angle of incidence of 1.58 from the normal. With this wave direction, Ly is 1200 m at Tp. The distance between the outer peaks of the symmetric mode shape no 6 with frequency ratio at Tp equal to 0.96, is 1100 m. It is thus natural that the bending moment response is now dominated by ampli®cation of mode 6. The maximum value away from the ends of the tunnel is now 1.35  106 kNm. According to our discussion the peak response may increase further for a wave direction of 1.758 from the normal, corresponding to Ly 1 1100 m. This is con®rmed with a maximum value of 1.55  106 kNm and similar distribution along the bridge as in Fig. 17. Short crested waves are de®ned with an energy spreading function C(y )=C cos8y where y is the angle relative to the mean wave direction (de®ned in the range 2p/2) and C is a constant determined according to energy preservation. With a mean wave direction normal to the vertical plane of the submerged tunnel, the maximum expected peak response away from the ®xed ends has dropped to 0.65  106 kNm compared with 1.85  106 kNm for long crested waves. Short crested waves with a mean direction 18 from the normal will give a maximum response equal to 0.57  106 kNm as compared with 0.50  106 kNm for long crested waves.

For a vertically curved submerged ¯oating tunnel, the present study has demonstrated high sensitivity of bending moment response to directionality for long crested waves. Response analyses for design purposes will have to consider a range of variation in the peak period Tp of the sea elevation spectrum. This is in order to ®nd the maxima related to close natural periods and to excitation from consecutive waves of varying non-normal incidence. It is also obvious that the actual curved form of the tunnel axis will in¯uence the response. For short crested waves there is a general reduction in response for this bridge and the directional sensitivity is not pronounced. In cases of counteracting directionality e€ects for long crested waves as for the 18 case above, the response may become higher for a short crested wave load model. An economical and still safe design will require reliable information about the actual wave spreading for the selected site of the submerged ¯oating tunnel. 4.1.3. E€ects of hydrodynamic mass, damping and geometric sti€ness for a vertically curved tunnel with tether mooring The vertically curved bridge moored with tension legs to anchors on the bottom depicted in Fig. 12 is chosen for evaluation of the e€ects of hydrodynamic mass, damping and geometric sti€ness, respectively. The analyses are carried out in the time domain with the general purpose non-linear FE program FENRIS [8]. This program has a highly modular structure which allows for combination of di€erent types of inte-

Table 3 Comparison of maximum sway moment computed at the midsection of the bridge Case

1

2

3

4

5

6

7

8

Sway moment 106 kNm

1.93

1.26

1.84

1.15

1.47

1.06

1.31

1.16

676

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

grated non-linear analysis, with particular capabilities for o€shore structures. As in the previous study of the vertically curved bridge in the frequency domain, the bridge is modeled by 50 equally sized 3-D beam elements. The tension legs are now modeled by cable elements, thus accounting for the true geometric sti€ness both in heave and sway. The waves are generated from a JONSWAP spec-

trum with signi®cant wave height, Hs=2.03 m, and with a peak period, Tp=5.5 s, using Fast Fourier Transforms of the wave spectrum. Long crested waves at a zero angle of incidence is assumed for all analyses in the time domain. The hydrodynamic loading and ¯uid interaction is modeled by way of the modi®ed Morison formula. The distributed force per unit length perpendicular to the length axis to the element, may then be written

Fig. 18. Variations in sway moment along the entire bridge obtained with the time domain analysis: (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (e) Case 5, and (f) Case 6.

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 18 (continued)

677

678

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 18 (continued)

fM ˆ

1 Cd rw D…un ÿ rÇ n † j un ÿ rÇ n j 2

…4:1†

‡ rw Aw uÇ n ‡ Cm rw Aw …uÇ n ÿ rÈ n † ‡ BÇr n Cd, Cm and B are the drag, the added mass and the wave damping coecients, respectively. D is the outer diameter of the cross-section Aw, and rw is the mass density of the water. un and uÇ n are the water particle velocity and acceleration component normal to the element axis, respectively, while Çrn and rÈ n are the struc-

tural velocity and acceleration component normal to the element axis, respectively. In contrast to the analyses carried out in the frequency domain, the frequency dependence of the hydrodynamic coecients Cd, Cm and B cannot be accounted for in the current version of FENRIS. In this study we therefore have kept these coecients constant throughout the entire simulation. The viscous e€ects are assumed to be small, hence the drag part is neglected by assuming Cd=0 for all elements.

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

An initial static analysis is carried out in one incremental step in which the gravity and the buoyancy loads are applied in order to establish the basis for the dynamic analyses. The dynamic analyses are performed by restarts from the static analyses. The time step

679

length is chosen to 0.25 s and the total simulation time is 600 s. In all analyses the irregular wave is gradually applied by prescribing the amplitude of the wave, and the ramp period is 10 s. Structural damping is accounted for by means of sti€ness and mass-pro-

Fig. 19. Time history plots for the midspan sway moment with the time domain analysis: (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (d) Case 5 and (f) Case 6.

680

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 19 (continued)

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 19 (continued)

681

682

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

p Fig. 20. Snapshot at time t=36.0 s for the 832-element model. The color fringes indicates the velocity magnitude v= u2x ‡ u2z m/s. The black horizontal line indicate the initial water surface z=0.

portional Rayleigh damping CT ˆ aM M ‡ aK KT

…4:2†

with damping coecients aM=aK=0.004. In order to study the di€erent e€ects, eight cases of dynamic analyses are considered (see Table 2). Fig. 18 displays the envelopes of the sway moment along the entire bridge for cases 1±6, see Table 2. The time history of the sway moment at the midsection of the bridge is shown in Fig. 19 for the same cases. Table 3 shows the corresponding maximum sway moment computed at the midsection of the bridge for all the considered cases. By comparing the results from the analyses it is seen that the modeling of damping is crucial. First we observe the big di€erence between the maximum sway moment in cases 1 and 2, 3 and 4, and, 5 and 6. For which the only di€erences is the modeling of wave damping, which is either neglected (cases 1, 3 and 5) or accounted for by picking the wave damping coecients from tables in ref. [11] for a frequency corresponding to the peak period Tp=5.5 s (cases 2, 4 and 6). Neglecting the wave-damping results in 39±60% higher values in the maximum sway moment at the midsection of the bridge. By comparing the di€erence in the sway moment in Table 3, between case 1 and 2, and, 5 and 6, it is seen that by increasing the structural damping the e€ect of neglecting wave damping is reduced from 53% to 39%. Furthermore by comparing the maximum sway moment obtained in case 1 and 5, and, 2 and 6, we see

that increasing the structural damping coecients by a factor of 2.5 from 0.004 to 0.01, results in a corresponding reduction of respectively 24% and 16%, depending on whether wave-damping is accounted for or not. We next consider the e€ect of hydrodynamic added mass. By comparing the results from Table 3 for case 1 and 3 (without wave-damping), and, 2 and 4 (with wave-damping), it is seen that picking the frequency dependent added mass coecients from tables for a frequency corresponding to the peak period or simply using the value Cm=1.0, for all elements, results in a 5% and 9% reduction, without and with wave-damping, respectively. The added mass coecients picked from the tables in [11] vary from Cm=0.43 for the elements connected to the endpoints to Cm=0.91 for the eight elements around the midpoint of the bridge. Case 3 and 4 (Cm=1) results in larger mass than case 1 and 2 (Cm=0.43±0.91), and a corresponding change in eigenmodes. This will cause a change from sti€ness to mass dominated dynamic ampli®cation, and thus explain the reduced ampli®cation factors obtained for the cases with Cm=1.0 compared to Cm=0.43±0.91. The e€ect of modeling the tension legs with non-linear cable elements instead of linear springs results in a reduction of the maximum sway moment of 4% (compare case 2 and 7 in Table 3). The net buoyancy of the tunnel is approximately 10% of the submerged weight for all cases but the last one for which it is reduced to 5%. This causes reduced tension in the tethers and thereby reduced geometric

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

683

Fig. 21. Force coecients for the cylinder.

sti€ness in the tension legs. By comparing the results from cases 2 and 8, for which the only di€erence is the net buoyancy, it is seen that the reduction in sway moments is reduced by approximately 10%.

In Fig. 21 the normalized drag and lift coecients obtained with SPECTRUM are plotted as functions of time. These are computed through, respectively. Cd ˆ

4.2. Numerical results for the Navier±Stokes ¯uid/ structure interaction problem The 2-D FSI-problem described in Section 3 is now simulated for a wave period of T=10.0 s. The wavenumber is then found from Eq. (3.5) as k=0.04027 mÿ1. The wave amplitude is increased linearly from Za=0 at time t=0 to Za=2.0 m at time t=T, from which it remains constant. We also apply a small constant current uc=0.5 m/s. A ®xed time step Dt=0.1 s is used in the simulation. Before the full FSI solution procedure is launched, an initial, pure CFD simulation is performed on a ®xed mesh (i.e. no solid region and no ALE stagger), starting from the initial conditions given in Section 3.1.1 and with a constant in¯ow velocity equal to uc=0.5 m/s. This simulation is continued until a nearly steady state is reached. The solution at the ®nal time step is then used as starting condition for the subsequent FSI-simulation, in which the water surface and the cylinder are free to move. An instantaneous picture of the ¯uid domain at time t=36.0 s (Step 360) is shown in Fig. 20. The gray tones indicate the absolute value of the ¯uid velocity.

Fx 1 r u2 DL 2 w c

and C1 ˆ

Fz ÿ rw gAw L 1 r u2 DL 2 w c

…4:3†

where Fx and Fz are respectively the horizontal and vertical components of the total force acting on the cylinder. Note that the buoyancy forces due to the hydrostatic pressure is subtracted from the vertical force, Fz, such that C1 is zero in still water. The displacement, velocity and acceleration of the cylinder are plotted in Fig. 22. For comparison, we have also performed a global dynamic analysis of the beam using simpli®ed models for the environmental loads. This analysis is performed with the structural FE code FENRIS where the wave load is calculated by means of the Morison formula. We see that in the vertical direction we get comparable results in FENRIS and in SPECTRUM for the velocity and acceleration. The amplitudes compare well also for the vertical displacement. However, it is also observed that a superimposed long-periodic mode, results in SPECTRUM. It seems to correspond with the natural period (45.5 s) of the structure, although the natural period is far above the wave period T=10.0 s. For the horizontal response there is a larger discrepancy between the FENRIS and SPECTRUM results.

684

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

Fig. 22. Comparison of global response parameters at the mid-span obtained with FENRIS and SPECTRUM. (a) Horizontal displacement (sway). (b) Vertical displacement (heave). (c) Vertical velocity. (d) Vertical acceleration.

As the structural properties of the FE model are equal in the two directions one should also expect response of similar magnitudes for the two directions, and this is the case for the FENRIS simulations. For SPECTRUM, however, the horizontal velocity and acceleration response is less than 50% of the vertical response. We have not yet been able to determine the cause for this di€erence.

5. Conclusions Methods for stochastic dynamic response analysis for submerged ¯oating tunnels subjected to wind driven surface waves are presented. The e€ect of directionality of the waves are demonstrated to be pronounced. Maximum levels of the response quantities may occur for non-normal incidence of the waves. This may be the case both for straight tunnels and for tunnels that are curved in the vertical or horizontal planes. Introduction of short crested

waves will reduce the directionality e€ects, and in general the response level will decrease with wave energy spreading. From the analyses carried out in the time domain it is demonstrated that the e€ect of modeling correct damping, both structural and hydrodynamic, is crucial for the vertically curved bridge moored with tension legs. The e€ect of modeling correct hydro dynamic added mass is much smaller than for the damping, but not negligible. It is demonstrated that the geometric sti€ness e€ect of the tension legs may be accounted for by means of linear springs. It is also noted that reducing the net buoyancy from 10% to 5% results in a corresponding reduction in the sway moment of approximately 10%. A possible reference solution for the ¯uid/structure interaction is formulated as a ®nite element implementation of the Navier±Stokes equations for the ¯uid region. Particular interest is related to structural components with typical 3-D ¯ow. Results for a section of a submerged tunnel in regular waves and current are

S. Remseth et al. / Computers and Structures 72 (1999) 659±685

presented. Comparison with a linear potential theory load e€ect analysis show good agreement for the vertical motion. There are still problems with a robust solution for the in-line motion of the bridge section, as predicted by the Navier±Stokes solver. Further developments are aiming at a practical interface between the interaction quantities computed by the FSI code and the dynamic response analysis codes applied for the global response analyses. Acknowledgements

[2] [3] [4] [5] [6]

The authors gratefully acknowledge the Norwegian Public Road Administration for their co-operation and support of research and development associated with submerged ¯oating tunnels. The study on structure/ ¯uid interaction using SPECTRUM was made possible through a grant from the Norwegian Research Council. The appreciation is extended to Vice President, Professor PaÊl G. Bergan of Det Norske Veritas for his co-operation and keen interest in this study.

[10]

References

[11]

[1] Borgman LE. Ocean wave simulation for engineering de-

[7] [8] [9]

685

sign. Journal of Waterways and Harbours Division, ASCE 1969;95(4):556±83. Shinozuka M. Monte Carlo solution of structural dynamics. Computers and Structures 1972;2:855. Hammersley F, Handscomb P. Monte Carlo methods. London: Methuen, 1964. Mathisen KM, Bergan PG. Large displacement analysis of submerged multibody systems. Engineering Computations 1992;9:609±34. Centric Engineering Systems Inc.: Palo Alto, California. SPECTRUM Solver (ver. 2.0) Command Reference and Theory Manual, May, 1993. Nomura T, Hughes TJR. An arbitrary Lagrangian± Eulerian ®nite element method for interaction of ¯uid and a rigid body. Computer Methods in Applied Mechanics and Engineering 1992;95:115±38. Faltinsen OM. Sea loads on ships and o€shore structures. New York: Cambridge University Press, 1990. Fyrileiv O. SESAM User's Manual, FENRISÐFinite element nonlinear integrated system. Hùvik, Norway, July: Det Norske Veritas Sesam A S, 1994. FEDAF, Dynamic analysis of ¯oating and submerged tubular bridges, User's Manual, SITNEF-report STF71 F88022, Trondheim, 1989. Hasselmann K, et al. Measurements of wind±wave growth and swell decay during the joint north sea wave project (JONSWAP). ErgaÈngzungsheft zur Deutschen Hydrographischen Zeitschrift, A No 12, 1973. Greenhow M. Hydrodynamic loading on horizontal cylinders. Unpublished project memo, Marintek, Trondheim, 1988.