1
Time domain simulations of wind- and wave-induced load effects on a three-span suspension bridge with two floating pylons 1)
Yuwang Xua 2), Ole Øisetha, Torgeir Moanb , c
a Department of Structural Engineering, Norwegian University of Science and Technology (NTNU)
bCentre for Ships and Ocean Structures (CeSOS), NTNU
cCentre for Autonomous Marine Operations and Systems (AMOS), NTNU
Abstract
The construction of a three-span suspension bridge with two floating pylons is currently being considered for crossing the 5-km-wide and 550-m-deep Bjørnafjorden in Norway. The bridge design represents a novel concept that requires a detailed dynamic analysis to improve the current understanding of its dynamic behavior. Geometric nonlinearities in the cables and mooring system and nonlinearities in the load models are of particular interest; in addition, the relative influence of the buffeting wind forces and the first- and second-order wave excitation forces were carefully studied.
The response predictions were obtained using state-of-the-art time domain methods.
Keywords: suspension bridge; floating pylon; wind load; wave load
1. Introduction
The Norwegian Public Roads Administration is administering a project—Ferry free coastal route E39—which aims to eliminate all ferries along the coastal highway E39 in Norway. One of the straits, Bjørnafjorden, is up to 5 km wide and 0.5 km deep, calling for a significant extension of current bridge technology. A three-span suspension bridge with two floating pylons, a combination of offshore and bridge technology, is a new concept for crossing wide and deep fjords [1, 2]. The bridge represents an entirely new design, requiring a detailed analysis of its dynamic behavior. Time domain methods are commonly applied when nonlinearities must be considered, as it is challenging to include such nonlinearities in a frequency domain analysis. Thus, assessing the influence of nonlinearities in the model is of particular interest to determine whether the calculations must be performed in the time domain or if the faster frequency domain methods are sufficient.
Modeling the motion-induced forces is a major challenge in the time domain simulations of the dynamic response, as they are dependent on the motion history. It is convenient to model the self- excited force in the time domain based on quasi-steady theory [3] and use coefficients from static wind tunnel experiments because the coefficients in quasi-steady theory are not dependent on the frequency;
however, it can be challenging to accurately model the self-excited forces using quasi-steady theory, which has resulted in a number of suggestions for improvements [4, 5]. The fluid memory effect can be considered by transfer functions in the frequency domain or by convolution integrals in the time domain. The self-excited forces for bridge decks are commonly modeled in the frequency domain using flutter derivatives, as proposed by Scanlan and coworkers [6]. The flutter derivatives represent an empirical generalization of the analytical expressions of the self-excited forces for airfoils (i.e., Theodorsen’s function) [7]. The Wagner function [8, 9] is the time domain counterpart to Theodorsen’s function, and this work has also been generalized for bridge applications. Time domain simulations of self-excited forces for bridge applications commonly start with an empirical expression for the transfer function in the frequency domain or the indicial functions in the time domain. The challenge is to fit the various models to the experimental data of the aerodynamic derivatives. The
1) Editor of this paper: John Niedzwecki, Texas A&M University, College Station, Texas, USA
2) Corresponding author.
E-mail address: [email protected] Published in Marine Structures
DOI: 10.1016/j.marstruc.2017.11.012
2 most common approach is perhaps the use of rational functions, also known as Roger’s approximation [10]. The majority of the aforementioned studies were either limited by simplified systems using still- air modes as generalized coordinates or detailed studies of the performance of the methodology considering a bridge deck section model. In bridge design, it is necessary to include self-excited forces in a finite element analysis of the entire bridge. Borri et al. [11] expressed the self-excited forces in the time domain by indicial functions and implemented the methodology in the finite element code FEMAS. Salvatori and Spinelli [12] also simulated the self-excited forces by convolution integrals using indicial functions, developed a finite element program capable of handling simplified bridge models, and analyzed the effects of structural nonlinearities and wind coherence. Chen et al. [13, 14]
used a state-space model to simulate the fluid memory effect, which can be more computationally efficient than solving the convolution integrals. Additionally, nonlinear effects were carefully studied in a flutter and buffeting analysis in the time domain. Øiseth et al. [15] also applied a state-space model of a simple beam with properties similar to a long-span bridge, and state variables were included as additional degrees of freedom in each node of a beam element.
The hydrodynamic radiation forces are similar to the aerodynamic self-excited forces in that they also depend on the motion history. The most convenient approach is to replace the frequency-dependent added mass and damping by constant coefficients, which are chosen at a dominating frequency, for instance, the peak frequency of the wave or the natural frequency of the structural system. However, this simple method cannot provide accurate results in the analysis of the structures’ transient response under a single frequency excitation or the steady-state response under multiple frequency excitations [16]. Cummins’ equation is widely used for time domain simulations of structures interacting with water to consider the frequency-dependent characteristics [17]. This equation is a vector integro- differential equation that involves convolution terms that account for the fluid memory effect and has been applied by many researchers [18, 19]. However, it is time consuming to solve the convolution integrals during a dynamic analysis [15, 20], and replacing the convolution integral with a state-space model is an attractive alternative. Taghipour et al. [20] verified that the same accuracy as obtained by solving Cummins equation directly can be obtained by replacing the convolution integral with a state- space model, and the calculations are approximately eight times faster. They also validated the methodology by comparing their results to experimental data for a flexible barge [21]. The state-space modal has also been used by many researchers in different areas [22, 23].
As outlined above, state space models are commonly used in the modeling of both hydrodynamic and aerodynamic self-excited forces. However, there are few studies on the performance of the methodology for modeling the dynamic behavior of structures subjected to both wind and wave actions. Thus, this paper provides a brief introduction to the state space modeling of self-excited aerodynamic and radiation forces; a description of the inclusion of state space models in a finite element model of a three-span suspension bridge with two floating pylons is provided. Studies on the dynamic behavior of the bridge consider first- and second-order wave excitations as well as the mean wind and linear and nonlinear buffeting forces. The influence of the nonlinear effects on the models is also carefully studied in order to present some general trends of what is important to include in the modelling of this novel bridge concept.
2. Dynamic response of a suspension bridge with floating
pylons
3 Fig. 1 Three-span suspension bridge with two floating pylons. Illustrated by Arne Jørgen Myhre, Statens
vegvesen
Fig. 1 shows a three-span suspension bridge with two floating pylons crossing Bjørnafjorden in Norway. The main cables are supported by two fixed pylons at each end of the bridge and two floating pylons in the middle of the bridge. The bottom part of the floating pylon is similar to tension leg platforms moored by four groups of tethers, providing a high stiffness in heave, pitch and roll. The water depth is 550 m and 450 m at the left and right floating pylons, respectively. To assess the dynamic behavior of the bridge, it is necessary to consider wind loading on the girder and pylons as well as wave loads on the floating pylons. The equation of motion can be written as
(1) (2 )
(t) ( ) ( ) ( ) (t) (t) (t) (t) (t)
Hydro Aero
s + s t + s+ h t = mean+ Buff + se + WA + WA± − Rad
F F
M u C u K K u F F F F F F
))))))))(
))))))( . (1)
Here, Ms, Cs and Ks symbolize the still-air mass, damping and stiffness matrix, respectively, and u represents the degrees of freedom of the finite element model. FAero represents the wind actions, which consist of a time-invariant component Fmean due to the mean wind velocity, a dynamic component FBuff due to turbulence in the wind field and self-excited forces Fse generated by the motion of the girder.
FHydro represents the wave actions, which consist of the radiation forces FRad induced by the motion of the submerged part of the pylons, the hydrostatic restoring stiffness Kh, and the first- and second- order wave excitation forces, FWA(1) and FWA(2 )± . Viscous drag damping forces on the submerged part of the floating pylons have not been considered since the aerodynamic damping is significant in the low frequency range where the hydrodynamic potential damping is close to zero. Possible excitations due to vortex shedding are out of the scope of this study, but are relevant for the design of the bridge girder, hangers and tethers. Possible vortex induced motions of the pylons might also be of relevance. We will describe the methods used to model the wind and wave forces in this chapter, with a particular focus on modeling the motion-induced forces in an efficient manner in the time domain to reduce the computational effort required.
2.1 Efficient modeling of motion-induced forces
2.1.1 Radiation forces
When a floating structure oscillates in waves or still water, it generates outgoing waves, resulting in oscillating fluid pressures on the surface of the body [24]. The integrated hydrodynamic pressures give rise to radiation forces, which are defined as follows for a single-frequency motion:
( ) ( )
Rad = h ω + h ω
F M u C u (2)
Here, Mh( )ω and Ch( )ω represents the frequency-dependent added mass and potential damping matrices, respectively. When the oscillation frequency tends to infinity, the damping converges to zero, whereas the added mass becomes constant and frequency independent, as shown in Fig. 2.
( ) ( ) ( )
( ) ( ) ( ) ( )
h h h
h h h h
ω ω
ω ω ω
= + ∞
= + ∞ =
M m M
C c C c (3)
Fig. 2 Typical frequency-dependent added mass and damping coefficients
4 The frequency-dependent added mass and damping are commonly obtained via potential theory, as discussed in the next chapters. The radiation force zRad, which only accounts for the frequency- dependent terms mh( )ω and ch( )ω , can be expressed in the frequency domain as
( ) ( ) ( )
Rad ω = ω ω
z u
G H G (4)
Here, H( )ω is the hydrodynamic transfer function, H( )ω =iωmh( )ω +ch( )ω , and ( )
Rad ω
Gz and Gu( )ω represent the Fourier transform of the radiation force and the velocity of the rigid body. By applying the inverse Fourier transform to Eq. (4), the forces can be expressed in the time domain as follows:
(t) (t ) ( )
Rad ∞ t t td
=
∫
−∞ −z h u (5)
Here, h(t) is the inverse Fourier transform of H( )ω . Eq. (5) consists of 6 6× =36 convolution integrals, and it is time consuming to solve the convolution integrals during a dynamic analysis, as noted by several authors [15, 20]. Replacing convolution integrals with state-space models is an efficient alternative. In this chapter, for brevity, we will only consider the convolution integral contributing to the radiation force in the i-direction due to motion in the j-direction, zij(Rad)(t) ∞h tij( t)uj( ) dt t
=
∫
−∞ − , as the other integrals are handled in the same manner.For a stable linear dynamic system, the relationship between the input uj(t) and the output zij(Rad)(t) can be characterized by an ordinary differential equation (ODE) with high-order derivatives [20]:
( ) 1 ( ) ( )
( )
1 1 1 0
1
1 1 1 0
(t) (t) (t)
(t)
(t) (t) (t)
(t)
n Rad n Rad Rad
ij ij ij Rad
n ij
n n
m m
j j j
m m m m j
d z d z dz
q q q z
dt dt dt
d u d u du
p p p p u
dt dt dt
−
− −
−
− −
+ + + +
= + + + +
(6)
The transfer function of the system can be obtained by taking the Laplace transform
( )
1
1 1 0
1
1 1 0
( ) ( ) ( )
( )
Rad
ij ij j
m m
m m
ij n n
n
Z s H s U s
p s p s p s p
H s
s q s q s q
−
−
− −
= ′
+ + + +
′ =
+ + + +
(7)
Here, Zij(Rad)( )s and U sj( ) are the Laplace transforms of zij(Rad)(t) and uj(t), and H sij′( ) is the transfer function. The relation between the transfer function in the Fourier and Laplace formations is
( ) i ( )
ij s ij
H s′ =ω=H ω . The degree of the polynomial in the denominator is often taken as one order higher than that of the polynomial in the numerator, n= +m 1. By choosing a proper order and applying the inverse Laplace transform, the relation between the output (radiation force) and input (velocity) in Eq.
(6) can be expressed as a state space model, which has been previously described in detail [20]:
{ }
(H) (H)
( ) (H)
(t) (t) (t)
, 1, 2, , 6
(t) (t)
c c j
Rad
ij c
u z i j
= +
∈
=
X D X E
Q X
(8)
where
0 0
1
1 1
2
2 2
(H) (H) (H)
3 2
1 1
(t) 0 0 0 0 0
(t) 1 0 0 0 0
0 1 0 0 0
(t) , ,
0
0 0 0 0
0 0 0 1 1
(t)
T
c c c
n
n n n
q p
X
q p
X
q p
p q
X q p
−
− −
−
−
−
= = = =
−
−
X D E Q
，
X, D(H)c , E(H)c and Q(H)c are different among each of the 36 state space models and can be determined by curve fitting the expression presented in Eq. (7) to the transfer functions defined by added mass and damping coefficients. Considering all degrees of freedom and the frequency-independent term, the radiation forces can be written as
5
( ) ( ) ( )
1 2 6
( ) 6 ( )
1
(H) (H)
( ) (H)
( ) (t) (t)
(t) (t) (t) (t)
(t) (t)
(t) (t) (t)
(t) (t)
Rad h Rad
Rad Rad Rad T
Rad
Rad Rad
i j ij
c c j
Rad
ij c
z z z
z z
u z
=
= ∞ +
=
=
= +
=
∑
F M u z
z
X D X E
Q X
(9)
2.1.2 Aerodynamic self-excited forces
Similar to the radiation forces, the aerodynamic self-excited forces q= qy qz qqT can be expressed as follows:
( ) ( )
ae K ae K
= +
q C u K u (10)
Fig. 3 Aerodynamic forces acting on the bridge section considered in this study
The positive direction of the displacement u and the forces q are displayed in Fig. 3. The aerodynamic damping matrix,Cae( )K , and the aerodynamic stiffness matrix, Kae( )K , contain 18 aerodynamic derivatives, Pn*, H*n and An*, n∈
{
1, 2,..., 6}
,which are cross-sectional properties that are functions of the reduced frequency of motion K=(Bω) /V .* * * * * *
1 5 2 4 6 3
* * * 2 2 * * *
5 1 2 6 4 3
* * 2 * * * 2 *
5 1 2 6 4 3
1 1
2 , 2
ae ae
P P BP P P BP
VKB H H BH V K H H BH
BA BA B A BA BA B A
ρ ρ
= =
C K
Here, V represents the mean wind velocity, ρ is the air density, B is the width of the girder, and ω refers to the oscillation frequency of the bridge deck. The aerodynamic derivatives are commonly determined by wind tunnel tests [25]. Eq. (10) is only valid for a single-frequency harmonic motion, but the model can be extended to any periodic or aperiodic motion through Fourier integral representation:
( ) ( ( ) ( )) ( )
( ) ( )
q ae ae u
u
ω iω ω ω ω
ω ω
= +
=
G C K G
F G (11)
Here, Gq( )ω and Gu( )ω are the Fourier transforms of self-excited forces and displacements, respectively. By applying the inverse Fourier transform to Eq. (11), the self-excited forces in the time domain can be expressed as a convolution integral as follows:
(t) ∞ (t t) ( )dt t
=
∫
−∞ −q f u
Here, f(t) and u(t) are the inverse Fourier transforms of F( )ω and Gu( )ω , respectively. The aerodynamic derivatives are only available at discrete reduced frequencies. A rational function [13-15, 26], as expressed in Eq. (12), is widely used to curve fit the scattered experimental data of aerodynamic derivatives and calculate the integration in the equation above.
3 2
1 2 3
1
1 /
( ) ( )
2 /
N l
l l
i B i B V
V V i B V d
ω ω
ω ρ
ω
−
= +
= + +
∑
+F a a a (12)
The unknowns a1, a2, al+3and dl in the rational function can be obtained by a least squares fit to the experimental data of the aerodynamic derivatives.
6 Similar to modeling radiation forces, the transfer function presented above can be expressed as a state space model. Integration of the distributed self-excited forces by applying the principle of virtual work and introducing the state space model yields the following expression for the nodal forces [15, 27]:
1 2
( ) ( )
( )
(t) (t) (t)
(t) (t) (t)
(t) (t)
se se
ae ae
c c
ae
se c
= + +
= +
=
F A u A u z X D X E u
z Q X
(13)
where
2 2
1 1 2 2
0 0
1 ( )
( ) 1 ( )
4 5
1 2 3
3
1 1
( ) ( ) ; ( ) ( ) ;
2 2
[ ] ;
; [ ] ; .
[ ]
L T L T
ae T
c
ae ae
c c N
T N N
V y y dy V B y y dy
V d
V d B
d
ρ ρ
−
−
= =
=
= = ⋅⋅⋅
= ⋅⋅⋅
∫ ∫
A N a N A N a N
I E I I I
D I Q A A A
X x x x I
Here, 3 2 3
0
1 ( ) ( )
2
L T
l+ = ρV
∫
y l+ y dyA N a N and l (t) d Vl t (d V Bl / )(t ) ( ) d B e
t t t
− −
= −
∫
−∞x u u . The matrix N(y) includes
the shape functions, and L refers to the length of the beam element.
2.1.3 Implementation in ABAQUS
The equation of motion for the combined structure and flow system is obtained when Eqs. (9) and (13) are introduced into Eq. (1):
(1) (2 )
2 1
( ) ( )
( ) ( )
( 1
( ( )) (t) ( ) ( ) ( ) ( ) (t) (t) (t) (t) (t)
(t) :
(t) (t) (t)
(t) (t)
(t) : (t)
s h s s h se Rad Mean Buff WA WA
se
ae ae
c c
se ae
c
Rad R
Rad
t t
Aerodynamic part
Hydrodynamic part z
+ ∞ + − + + − − + = + + + ±
= +
=
=
M M u C A u K K A u z z F F F F
z X D X E u
z Q X
z z
) ( ) ( )
2 6
( ) 6 ( )
1
(H) (H)
( ) (H)
(t) (t) (t)
(t) (t)
(t) (t) (t)
(t) (t)
ad Rad Rad T
Rad Rad
i j ij
c c j
Rad
ij c
z z
z z
u z
=
=
= +
=
∑
X D X E
Q X
(14)
The motion-induced hydrodynamic and aerodynamic forces are modeled using similar state space models.
(t) (t) (t)
(t) (t)
c c
c
= +
=
X D X E u z Q X
(15)
As shown above, the velocity of the structure is the input of the model, whereas a vector of the motion-induced forces is the output of the model. X represents the aerodynamic or hydrodynamic state variables. These variables must be included when solving the equation of motion for the combined structure and flow system. The starting point is to transform the continuous first-order linear inhomogeneous differential equation presented in Eq. (15) to a discrete form [28].
1 [(1 ) 1]
k k k k k
k k
γ γ
+ = + + − + +
=
X DX Eu G u u
z QX
(16)
Here, D=eDc∆t, E=(D I D− ) c−1Ec, G=Dc−1(E E− ∆c t) and Q=Qc.
Algorithm 1 illustrates how the state variables are updated at each time step in the dynamic analysis by defining a user element. The elements are modeled on top of the nodes of ordinary beam elements such that it is not necessary to involve the mass, damping and stiffness terms related to the structure. A first-order hold assumption is introduced for the discretization of velocity; otherwise, G will be zero.
Under the first-order hold assumption, a much larger time step is allowed in the numerical integration of the equation of motion than can be used based on the zero-order hold assumption [27]. However,
7 the assumption makes the state-space model more complicated because the state variables are dependent on acceleration in the current time step. Therefore, it is necessary to reorganize Eq. (16) and introduce a new variableXk =DXk−1+Euk−1+G(1−γ)uk−1 such that Xk depends only on the motion in the previous time step and the term G uγk can be added directly into the mass matrix of the system.
h( )∞
M , A2 and A1 are constant matrices and can be conveniently modeled by defining them as the mass, damping and stiffness of the user element, respectively.
Algorithm 1. Simulation of the self-excited force by the user element in ABAQUS Step (1): ABAQUS supplies current estimates of uk, uk and uk
Step (2): Load Xk−1, uk−1 and uk−1 from the last time step (the initial value of X :X0=0) and calculate the state vector Xk−1
Xk−1=Xk−1+G uγk−1
Step (3): Calculate the residual force and Jacobian matrix Rk =Museruk+Cuseruk+Kuseruk+QXk
user 12 user user
t t
γ
β β
= + +
∆ ∆
J M C K
where Muser, Cuser, and Kuser are the mass, damping and stiffness, respectively, of the user element:
Muser= −QGγ ,Cuser = −A2,Kuser = −A1 for the aerodynamic user element Muser=Mh( )∞ −QGγ ,Cuser=0, Kuser =0 for the hydrodynamic user element Xk depends only on the motion in step k-1:
Xk =DXk−1+Euk−1+G(1−γ)uk−1
β and γ are the two parameters in Newmark’s method
Step (4): ABAQUS calculates new values for uk, uk and uk and performs equilibrium iterations until the results converge.
Step (5): Save Xk, uk and uk in ABAQUS as solution-dependent variables Step (6): Return to (1) for new time step or stop
2.2 Buffeting forces
2.2.1 Simulation of the wind velocities
The wind field is simulated using Monte Carlo analysis [29, 30]. The density of the air is assumed to beρ=1.25kg/m3, and the cross-spectral densities of the horizontal along-wind velocity u and vertical components w at points i and j are assumed to be given as follows:
2 2
5/3 5/3
5/3
7/3
40.58 3.18 ( ) ( )
40.58
( , , , ) exp( )
(1 9.74 / ) (1 9.74 / )
( , ) 0.82 exp( )
(1 0.79 / )
( , ) 2.23 exp( )
(1 1.67 / )
j j i i
uu
i i j j i j
ww
uw
V z x z
S x z z V z
z V z V V V
Vz x
S x
z V V
Vz x
S x
z V V
κ ω
ω κ
ω ω
κ ω
ω ω
κ ω
ω ω
+
+
+
∆ + ∆
∆ ∆ = −
+ + +
∆ = −∆
+
∆ = −∆
+
(17)
where κ is the roughness coefficient at the site, which is assumed to be 0.0031; zi and zj are the heights above the ground at the two points; and ∆x and ∆z are the lateral and vertical distances, respectively, between the two points considered.
The vertical curvature of the girder is neglected by using the height at the middle of the bridge as the reference for the wind load of the girder. It is further assumed that only the along-wind turbulence component u provides wind loading on the pylons and that the vertical mean wind velocity profile is given by
8 ( )
ref ref
V h h
V h
α
=
(18)
Here, Vref and V h( ) are the mean wind velocity at the reference point and at height h, respectively, and α is related to the surface roughness. The wind velocities on the girder and pylons are simulated independently, neglecting a possible correlation. This assumption is clearly an approximation in the areas near the connection between the pylons and girder; however, previous results [31] indicate that the approximation is well within the overall uncertainty of the model considering that wind tunnel tests of the pylons are also currently unavailable.
2.2.2 Linear and nonlinear wind forces
The mean and buffeting forces due to the mean and turbulent wind are calculated using quasi-steady theory [32] when the aerodynamic admittance is neglected. Assuming that the fluctuating flow components, u x t( , ) and w x t( , ), and the structural velocity are small compared to the mean wind velocity V, the linearized wind-induced forces can be defined as
2 ( ) 2( ) ( )
( , )
(t) 2 ( )
( , )
2 2
2
D D D L
mean Buff L L L D
M M M
D B C D B C D B C C
u x t
V B VB
C C C D B C
w x t
BC BC BC
ρ ρ ′ ′ − ′
+ = + + ′
F F (19)
Here, D is the height of the girder. CD, CL and CM are the mean values of the drag, lift and torsional moment force coefficients, respectively, and C′D, CL′ and CM′ are their derivatives with respect to the angle of attack.
At the peaks of the turbulent wind velocities, the higher-order terms can significantly contribute to the wind loading and can thus significantly contribute to the load effects in certain situations. The buffeting forces are modeled using the following expression to investigate the influence of the nonlinear terms.
2
2
2 2 2
( )
cos sin 0
(t) 1 sin cos 0 ( )
2 0 0 1 ( )
( ( , t)) ( , t) ( , t)
arctan( )
( , t)
D D
mean Buff rel L L
M M
rel
D C C
V B C C
B C C
V V u x w x
w x V u x
β β β
ρ β β β
β
β
+ ′
−
′
+ = +
+ ′
= + +
≈ +
F F
(20)
2.3 Wave excitation forces
The hydrodynamic actions consist of hydrostatic forces, radiation forces and wave excitation forces [33]. The hydrostatic forces are frequency independent and modeled as springs in the finite element model of the structure, whereas the radiation forces are modeled using state space models, as presented in chapter 2.1. The wave excitation forces related to the fluctuating pressure on the submerged surface of the pylons are caused by incident waves assuming that the structure is fixed. As noted above, the fluid-structure interaction effects are commonly modeled using potential theory; in this paper, we will consider second-order difference and sum frequency forces as well as first-order forces. Higher-order models exist but are outside the scope of this paper.
2.3.1 First-order wave excitation forces
The first-order forces are always larger than the second-order forces and are defined by a transfer function that can be obtained by potential theory. The amplitude of the first-order wave forces is proportional to the wave amplitude. In irregular waves, the first-order forces can be obtained by summing the contributions from all of the frequency components and directions using the following well-known expression [34]:
9
[ ] { }
(1) (1)
(1) 1
(1)
( , , ) ( , ) cos ( cos( ) cos( )) , 1, 2, , 6
2 ( , )
Im( ( , ))
tan Re( ( , ))
N M
WA n m nm n m m n mn mn
n m
nm n m
n m
mn
n m
x y t k x y t i
Sη
ω q η q q ω e φ
η ω q ω q
φ ω q
ω q
−
= + − + − ∈
= ∆ ∆
=
∑∑
F T
T T
(21)
Here, T(1)( , )ω q is the first-order transfer function and represents the force generated by a unit amplitude regular wave with frequency ω propagating from the directionq . The wave amplitude and wave number are denoted by η and k, respectively; for small-amplitude waves, the wave number is related to the frequency and water depth h by the dispersion relation, ω =2 gktanh(kh). emn∈
{
022π}
is a uniformly distributed random phase angle. Sη( , )ω q represents the wave spectrum and is a function of both the wave frequency and direction.( , )= ( ) ( , )
Sη ω q S ω Dω q (22)
Here, S( )ω is the unidirectional wave spectral density and D( , )ω q is the directional distribution. The directional distribution for locally wind-generated sea is commonly approximated as frequency independent. Little information regarding the wave conditions at the site is currently available, making the one parameter Pierson-Moskowitz spectrum [35] a tempting choice. Meanwhile, cos-2s distribution [36] is applied to represent the wave direction distribution in this study.
2 5
4 2
2
( ) 0.0081 exp( 3.11 )
(s 1)
( ) cos
2 (s 1 / 2) 2
s s
S g
H D
ω ω
ω q q
π
= − −
Γ +
= Γ +
(23)
where Hs is significant wave height and s is a parameter which determines the concentration of the direction distribution of the wave energy.
2.3.2 Difference frequency forces
The natural frequencies of the low-order modes of the bridge are well below the energy content in the wave spectra such that these modes will not be excited by first-order excitation forces. However, the difference frequency forces can cause wave excitations in this frequency range. For short-crested waves, the different frequency forces can be written as follows [37-39]:
( ) t ( )
(2 ) (2 )
1 1 1 1
Re ( , , , ) k j kh jl
N N N N
i i
WA jl kh j k l h
l h j k
e ω ω e e e
η η ω ω q q − −
− −
= = = =
=
∑∑∑∑
F T (24)
where T(2 )− contains the in-phase and out-of-phase components of the full quadratic transfer function, Tic and Tis, respectively, and represents the forces induced by the interaction of a unit-amplitude wave associated with frequency ωj and direction ql and a unit-amplitude wave associated with frequency ωkand direction qh. If the direction interaction effects are ignored, i.e., the terms where
l h
q q≠ , the difference frequency forces can be written as
( ) t ( )
(2 ) (2 )
1 1 1
Re ( , , ) k j kl jl
N N N
i i
WA jl kl j k l
l j k
e ω ω e e e
η η ω ω q − −
− −
= = =
=
∑∑∑
F T (25)
This will reduce the computational effort significantly, as fewer quadratic transfer functions are required. The mean drift force is included in the equation above and can be found by setting ωk =ωj:
1 1
2 ( , ) ( , , , )
N N
ic
WA j l j j l l
l j
Sη ω q ω ω q q ω q
= =
=
∑∑
∆ ∆F T (26)
2.3.3 Sum frequency forces
The sum frequency forces introduce excitation at higher frequencies compared to the first-order excitation forces and may be important for the heave, pitch and roll motions of the pylons. If the directional interactions are ignored, the time series of summed frequency forces can be written as follows [39].
( ) t ( )
(2 ) (2 )
1 1 1
Re ( , , ) k j kl jl
N N N
i i
WA jl kl j k l
l j k
e ω ω e e e
η η ω ω q + +
+ +
= = =
=
∑∑∑
F T (27)
10
3. Numerical simulations
A comprehensive finite element model of the bridge, as displayed in Fig. 4, is used in the dynamic analysis. Each of the main girder spans has a length of 1385 m and the pylons are 198 m high and have a draft of 65 m. The bending stiffness of the girder for horizontal, vertical and torsional deflections are 3.22 10⋅ 12Nm2, 2.04 10⋅ 11Nm2, and 1.99 10⋅ 11Nm2 respectively. The main cable and hanger have a cross section area of 0.333 and 0.032 m2, respectively. More details about the finite element model of the bridge are presented in [40]. The girder, main cable, tethers, hangers and pylons are modeled by beam elements. The aerodynamic self-excited and hydrodynamic radiation forces are simulated by implementing a user element in ABAQUS. The user element is developed as a one-node element and is included in the nodes of the girder and a reference point on each of the floating pylons; the user elements are illustrated by the red markers in Fig. 4. The added mass, potential damping and first- and second-order transfer functions for the wave excitation forces are obtained using potential theory. Fig.
5 shows the hydrodynamic panel model of the submerged part of the pylon and the free water surface.
The panel model of the free water surface is only required when calculating the full quadratic transfer functions based on the near-field formulation [38]. The side span on the right hand side of the Artist’s view in Fig. 1 has not been included in the modeling.
Fig. 4 Finite element model of the three-span suspension bridge with two floating pylons
Fig. 5 One quarter of the panel model of the submerged part of the pylon and the water surface in WADAM Pylon 1#
Pylon 2#
Pylon 4# Pylon 3#
x y z
1# 2#
3#
4# Tethers
Pylon 2#
roll pitch
heave
yaw
x
y
z
11
3.1 Estimation of the parameters in the state space models
3.1.1 Potential added mass and damping
The potential added mass and damping at discrete frequencies are obtained by linear potential theory using the software WADAM [41]. The numerical results for selected components of the added mass and potential damping are displayed in Fig. 6. The rational function presented in Eq.(6) is used to curve fit the data. The unknowns θ=[pn−1,2,p q0, n−1,2,q0]T can be obtained by regression in the frequency domain:
(H) 2
2
(i , ) ( ) (i , ) arg min
(i , )
j l l ij l
l j l
U H Z
U
ω ω ω
ω
=
∑
θ − θθ θ
(28)
This method has been implemented in a MATLAB toolbox, FDI_Toolbox_v1.2, by Perez and Fossen [42], which is applied in this study. A high-order polynomial of the transfer functions must be applied, as there are many fluctuations in the numerical results for the added mass and damping. Fig. 6 shows that the models fit the data well.
Fig. 6 Selected added mass and damping coefficients. The order of the transfer functions n is given in the figures.
3.1.2 Aerodynamic derivatives
The self-excited forces are modeled by fitting rational functions to the transfer function defined in terms of the aerodynamic derivatives. The accuracy of the numerical simulations depends on the quality of the curve fit of the expressions presented in Eq.(12) to the experimental data of the aerodynamic derivatives. The self-excited forces must be captured precisely throughout the reduced frequency range, particularly in the range corresponding to the natural frequencies of the system.
Experimental aerodynamic derivative data [43] from forced vibrations tests of the cross section of Hardanger Bridge are represented by circles and triangles in Fig. 7, and the solid lines are curve fits defined by the rational functions introduced in Eq. (12). As shown in Fig. 7, the curve fits match well with the experimental data. We have applied the shape of the girder of the Hardanger Bridge in this study since wind tunnel tests are available for this section.
3.2 Simulation of the wind and wave loads
Eqs. (17)-(27) are used to generate the wind and wave forces acting on the pylons and girder. The surface roughness related parameter α in Eq.(18) is assumed to be 0.12 [31]. The fluctuating wind velocities have been simulated at 151 points along the girder and 40 points along the pylons, with cut- off frequencies of ωu =60 rad/s and ∆ω=0.001 rad/s, respectively. Rather small difference of the response statistics is observed when applying more points along the girder, which infers that 151 points can represent the wind field very well. The following force coefficients of the girder were
12 Fig. 7 Curve fits of the 18 aerodynamic derivatives obtained by forced vibration tests
Fig. 8 (a) Wind force in the y-direction over one span length (1,385 m). (b), (c) and (d) show the first-order, second-order difference and sum frequency wave excitation forces in the y-direction on one floating pylon, respectively.
obtained from wind tunnel experiments [44]: CD=0.70, CD′=0, CL=-0.25, C′L=2.4, CM=0.01 and C′M
=0.74. The wave amplitude of the first- and second-order wave excitation force transfer functions,
(1)
Twave and Twave(2 )± , are calculated by WADAM. The time series of the buffeting wind force and the first- and second-order wave excitation forces in the y-direction are shown in Fig. 8; here, a mean wind velocity of 35 m/s and significant wave height of 4.88 m are applied. The average estimated spectral
13 Fig. 9 Auto spectral density of the wind and wave forces presented in Fig. 8. Dashed lines are the target spectral densities, whereas each solid line represents the average spectral density for 10 realizations.
density of 10 realizations is compared to the target spectral density in Fig. 9 to verify that the time domain simulations correspond to the target spectral densities.
The magnitude of the wind force and second-order difference force is very small compared to the others, but they are dominant in the frequency range of 0 to 0.3 rad/s. The natural frequencies of the first ten modes and first three torsional modes considering the mean wind force at V=35 m/s and the added mass as the frequency goes to infinity are listed in Table 1. The first three natural modes of the structure have frequencies in this range, making the forces important for the dynamic behavior of the system. The first-order wave excitation forces are dominant in the frequency range of 0.3 to 1.2 rad/s, which may excite several higher modes of the girder. The second-order sum force is in the range of 1.7 to 3 rad/s, which can be important for modes where the pylons move in heave, pitch and roll.
Table 1 Natural frequency for the first ten modes and the first three torsional modes. H, V and T refer to horizontal, vertical and torsional vibration modes, respectively.
Mode no. Natural frequency (rad/s) Type
1 0.0724 H
2 0.0968 H
3 0.186 V
4 0.307 H
5 0.317 H
6 0.359 H
7 0.464 H
8 0.518 H&V
9 0.563 V
10 0.643 H
51 1.845 T
52 1.868 T
54 2.001 T
3.3 Computational efficiency and validation
The time domain methodology outlined in this paper should provide the same response as a multimode frequency domain analysis if the system is linear and a sufficient number of modes are included in the model. The time domain simulations should be verified by comparing the response statistics with results from a frequency domain analysis. Multimode theory is briefly introduced in the appendix.
