# Journal © 1997ElsevierScienceB.V.All rights reserved. PII S0 167-61

Journal of Wind Engineering ~ 1 ~

ELSEVIER and Industrial Aerodynamics67&68(1997)239-252

Numerical simulation of flow around a box girder

of a long span suspension bridge

Shinichi Kuroda 1

lshikawajima-Harima Hea~y Industries Co., Ltd., 19-10, Mohri 1-Chome, Koto-Ku, To~,o 135, Japan

Abstract

A numerical simulation of a high Reynolds number flow around a fixed section model with

a shallow hexagonal cross-section is presented. The two-dimensional incompressible

Navier-Stokes equations are used as the governing equations. The numerical algorithm is

based on the method of pseudocompressibility and uses an implicit upwind scheme. The overall

characteristics of the measured static force coefficients were well captured by the present

computation.

Keywords. Numerical simulation; Navier-Stokes equations; Method of pseudocompressibility;

Upwind scheme; Long span bridge; Box girder; Great Belt East Bridge

1. Introduction

At the present, wind tunnel testing is actually the only way to assess the aerodynam-

ic performance and aeroelastic stability of long span bridges. The full bridge model

test may be the final phase of the wind resistant design. Section model tests are also

used to improve the aerodyne.mic/aeroelastic performance of the bridge components

such as bridge decks or tower shafts. In recent years, specific analytical methods such

as the direct flutter FEM analysis 1 have been developed and used to complement

the full bridge model test. However, in these analyses, static/dynamic section model

tests are mandatory to obtain the wind force data, which are the input for the analyses.

Even a section model test is generally very time-consuming and expensive. With the

recent rapid improvement of computers, Computational Fluid Dynamics (CFD) has

made a remarkable progress. If the section model test can be replaced by the CFD, it

will be significant for bridge engineers.

The objective of the present research is to develop a numerical method to simulate

the static/dynamic wind tunnel test for section models of long span bridges. The

1E-mail: [email protected]

0167-6105/97/$17.00 © 1997ElsevierScienceB.V.All rights reserved.

PII S0 167-61 05 (97)00076-7

JOURNAL OF

240

S. Kuroda/’J. Wind Eng. Ind. Aerodyn. 67&68 (1997) 239 252

H4.1

o

75

75 225

O- 375 (30.0)

Fig. 1. Section model 2.

27.0 m

Fig. 2. Box girder cross section of the Great Belt East Bridge 3.

simulation of the static wind tunnel test is the first step of this CFD research. In this

paper, a numerical simulation of flow around a fixed section model with a shallow

hexagonal cross-section 2 (Fig. 1) is presented. The cross-sectional configuration

was a candidate for the girder section of the Great Belt East Bridge 2,3 (Fig. 2).

Extensive section model tests, as well as full bridge model tests, for the Great Belt East

Bridge were conducted at the Danish Maritime Institute (DMI) and the static force

coefficients for the section model of Fig. 1 are found in Ref. 2. The computed static

force coefficients were compared with this data.

2. Numerical procedure

The governing equations are the incompressible two-dimensional Navier-Stokes

equations described in a generalized coordinate system. The numerical scheme em-

ploys a finite-difference procedure and the equations are solved using the method of

pseudocompressibility I-4.

In this paper, only the stationary grid was used, however, when considering the

translation/deformation of the coordinate system itself, the generalized coordinates

are as follows:

= ~(x, y, t), (la)

= ~(x,y, t). (lb)

H – 54.7 (4.38)

J.7

37.5

S. Kuroda/~ Wind Eng. Ind. Aerodyn. 67&68 (1997) 239 252 241

Using these, the governing equations are written in conservation form as follows:

– + ~

(?t- 04

= 0,

~ (f-f)

1~gxP+ uU + {,u1

= -J Lg,P + vU + ~,v J’

/=1! + uv + .,u-1

I

krlyp + t,V + rhVJ”

(2a)

= – j*” (2b)

,2c,

(2dl

(2e)

a=j ,

1

Eq. (2a) represents the conservation of mass and Eq. (2b) the conservation of x and

y components of the momentum, respectively. J is the Jacobian of the coordinate

transformation, and U and V contravariant components of the velocity vector. ~x, etc.

are the metrics of the transformation.

j- 1 __ x~y, – ycx,, (3a)

U = ~xu + ~.v, (3b)

V = ~lxu + ~,v, (3c)

¢~- dx’ etc. (3d)

The viscous fluxes are written as follows with the kinematic viscosity as v:

~ = j (~2 + ~2)v¢ + (~rL, + ~yrL,.)v,’ (4a)

(4b)

In the case of a stationary grid, the generalized coordinates are as follows:

= ~(x, y), (%)

~l = it(x, y). (5b)

242 S. Kuroda/J. Wind Eng. lnd, Aerodyn. 67&68 (1997) 239-252

The time derivatives of 4 and q are zero (4, = q, = 0). The inviscid fluxes become

1 ~_xP4-uU1,

= J Ld,,p+ vUJ (6a)

l lq-~p+uV1 (6b)

f =- J Lrlyp + vV ”

Other terms are not changed because the time derivatives of ~ and r/do not appear

in them. In this paper, Eqs. (6a) and (6b) were used for the inviscid fluxes.

Using the pseudocompressibility method 4, subiterations are required to satisfy

the continuity equation at each physical time step. For evaluating the convective

terms in the right-hand side, the upwind differencing scheme based on the flux-

difference splitting 4 was used. The upwind-biased scheme of third order or fifth

order of Refs. 4,5 can be used. For all cases presented in this paper, the fifth order

was used. For evaluating the viscous terms, the second-order central differencing was

used. To solve the algebraic equation, the implicit line-relaxation scheme was used.

Implicit boundary conditions are used at all of the boundaries. At the wall boundary,

the no-slip condition was adopted for the velocity, and for the pressure, such a condi-

tion was imposed that the component perpendicular to the wall surface of the pressure

gradient becomes zero. For the inflow and outflow conditions, characteristic bound-

ary conditions 4 were used.

3. Conditions of computation

In Ref. 2, no description can be found of the wind speed or Reynolds number at

which the static force measurements were conducted. However, the following descrip-

tion can be found in Ref. 3. Steady state wind load coefficients were measured on

1 : 80 and 1 : 300 models for the selected section H9.1 which was identical to the

section of Fig. 1 (H4.1) except that the deck width of H4.1 was 1 m narrower than

H9.1 (30 m versus 31 m for H9.1). Satisfactory agreement was demonstrated between

steady state wind loads obtained from models of different size, indicating that

Reynolds number effects are of minor importance in the case of H9.1.

In the present computations, a Reynolds number of 3 x 105 (based on the girder

width B) was used, assuming a wind speed of around 12 m/s. The experimental static

force measurements were conducted for angles of attack ranging from – 10° to + 10

2. Thus, the computations were carried out at six angles of attack ~ of – lff, – 5′,

0°, Y’, 8° and 10°. The measurements were conducted in turbulent flow with a longitu-

dinal turbulence intensity of 7.5% 2. However, in the present computations, a uni-

form smooth flow was assumed as the incoming flow and no turbulence model was

used. The turbulent flow simulation together with the three-dimensional extensiton

will be presented in a future work.

The O-type grid (221 x 101) was generated around the girder surface using the

hyperbolic grid generation scheme 6. A close view of this grid is shown in Fig. 3.

S. Kuroda/J. WindEng. Ind. Aerodyn. 67&68 (1997) 239 252 243

Fig. 3. Grid used in the present computations.

Identical grids were used for all angles of attack. The normal spacing next to the wall

was 0.0002 girder widths. The outer boundary of the computational region was

placed about 6 times the girder width away from the girder in all directions. The

presence of side railings and central crash barrier was not considered in the present

computations. The inclusion of such equipments in the computation will be also

presented in a future work.

A non-dimensional time step At (= AT U/B) of 0.01 was used. The sweep for the

line-relaxation scheme was 1 time per subiteration step each for ~ and t/directions. As

the convergence criteria for the iterative calculations, ipm+1 _pm I < 10- 3was used for
the pressure, and for the momentum equation, the residual of x or y components was
not more than 10- 4. In the above, the index m indicates the subiteration count at each
time step.
4. Computed flow field
The computations were started impulsively from the uniform flee-stream condi-
tions and were continued until the non-dimensional time t ( = T U/B) equal to 40.
Only for 2 = 10°, the simulation was continued until t = 80, because the period
of flow field variation (although not exactly periodic) was long compared with those
for other angles of attack. Fig. 4 shows the computed instantaneous streamlines at
t = 40.
2 4 4
S. Kuroda/J. Wind Eng. Ind. Aerodyn. 67&68 (1997) 239 252
= _10 °
4.1. Meanflow pattern
~ = 5 °
a=O °
a=lO °
=-5° ~=8°
Fig. 4. Instantaneous stream lines.
The mean flow patterns are shown in Fig. 5. The averaging procedure was taken in
the range of t = 20~40 ( for c~= 10°, the range of t = 40-80).
c~= 5~',8°, 10°: On the lower surface, the flow is attached to the inclined web surface
on the windward side from the stagnation point to the corner. The flow shows a small
separation at the corner and reattaches the lower flange soon again and runs to the
corner on the leeward side. The small separation region on the windward side does
not vary in size in the range of these positive angles of attack.
fjJ
J
S. Kuroda/a~ Wind Eng. Ind. Aerodyn. 67&68 (1997) 239-252 245
ot = 0 °
J
j/f
a = 10 °
o~ = - 5 °
oz = 8 °
oz=-10°
J
Fig. 5, Time-averaged stream lines.
a =5°
On the upper surface, the separation occurs at the leading edge and the flow
reattaches on the surface of the upper flange. The higher the angle of attack, the larger
the size of the separation bubble originated from the leading edge.
= 0°: The separation does not occur at the leading edge and the flow attaches on
the upper and lower inclined web surfaces on the windward side. Flow separation
occurs at the corners. The size of the separation bubble on the bottom plate is larger
than that at the positive angles of attack (5°, 8°, 10°).
= - 5, - 10°: At both angles, the flow fundamentally attaches on the upper
surface. On the lower inclined web surface, the leading edge separation is small at
246 S. Kuroda/J. Wind Eng. lnd. Aerodyn. 67&68 (1997) 239 252
= - 5°, however, at ~ = - l0 °, large separation occurs at the leading edge and the
flow reattaches just at the corner. As for the flow over the bottom plate, the flow does
not reattach on that surface at cx= - 5~ but does at ~ = - 10°. The size of the
separation bubble on the bottom plate is almost the same between the cases of
=- 10°and~=0°.
4.2. Mean surface pressure
The mean surface pressure distributions are shown in Fig. 6. The pressure coeffic-
ient Cp in Fig. 6 is normally defined as
2p
Cr -- pU 2" (7)
Upper surface: At e = 5°, 8 and 10°, the pressure on the upper surface becomes
negative all over the surface. On the whole, the higher the angle of attack in the
positive direction, the lower the surface pressure. However, the value of the negative
peak at ~ = 8° is lower than that at ~ --- 10°. The peaks of the negative pressure are
located around the corner on the windward side.
At ~ = 0°, the pressure becomes negative from the middle on the inclined web
surface, and shows an almost constant (negative) value on the first half of the
separation bubble. On the last part of the separation bubble, the pressure rises toward
the reattachment point. Downstream from the reattachment point, the pressure is
almost constant.
At ~ = - 5, - 10°, the pressure on the upper surface is positive nearly all over the
surface except in the vicinity of the inclined web surface on the leeward side. At both
angles, the pressure displays almost the same distribution along the upper surface,
however, somewhat higher at c~= - 10° on the whole.
Lower surface: At c~-- 5, 8 and 10°, the pressure on the lower surface displays similar
distributions. With the acceleration of the flow along the inclined web surface, the
pressure decreases monotonically from the stagnation point. The pressure is negative
inside the separation bubble on the bottom plate and rises toward the reattachment
point. Downstream from the reattachment point, the pressure decreases slightly in
a steady fashion.
At e =0 °, the pressure distribution on the lower surface displays some similarity to
that at the positive angles of attack. However, on the whole, the value of pressure at
= @ is lower than that at the positive angles of attack, At c~= 0°, the pressure
becomes negative nearly all over the surface except in the vicinity of the leading
edge.
At ~ =- 5 and -10 °, the pressure is negative all over the surface. On the
inclined web surface on the windward side, the pressure at ~ = - 10° shows a far
lower value than that at c~= -5 °. However, on the last half of the bottom
plate, the pressure at ~ = -5 ° shows a negative peak and is lower than that at
~= -- 10 °.
Q.
(9
•
-° ..f-o..~_.
Q,.
(.b
-0.5oj
05
1
1.5
x/B
Fig. 6. Time-averaged pressure coefficient distribution.
S. Kuroda/J. Wind Eng. Ind. Aerodyn. 67&68 (1997) 239-252 247
Upper Surface
-2.5
-1.5~ -
f
i~-
--~= 10°
--~.~= 8 °
~..~= 0°
...... ~, =-I0"
-1
-0.5 .......~--
I
....
........................ ........
0.5 '/'L"
!#I
1.5 I
0 0.2 0,4 0.6 0.8 1
-2.5
-2
15
....
2
-I ....
....
x/B
Lower Surface
i
0 02 04 0.6 0.8 1
i'
~a=lO °
..... ~= 5u
..... ~,= .5°
..... a,=,10°
248
S. Kuroda/J. Wind Eng, Ind. Aerodyn. 67&68 (1997) 239 252
T5
~J
(~ 0.5
o
io, I
I
.15d
15
I
0.5
u.
i -o.5
=8 _,
-1.5
-2
5 10 15 20 25 30 35 40
TiME
0 10 ,20 ~ID 40 50 80 70 BO
TIME
~2
-os -os
g,
2
ls
2i-- --r, .....
15
i
...................
' .........
' ....
' ....
~- 0°
~=10"
a-.5
°
a'= 8"
0 5 10 15 20 25 t0 35 40
TIME
o.=.10 °
5 to 15 1~10 25 3o :15 40
TIME
S 10 tS 20 25 30 35 40
TIME
$ 10 15 20 25 30 35 40
TIME
--col
3
o
f~
ol
8 ~2
Fig. 7. Time histories of drag, lift and pitching moment.
C
L
CM
2
0.5
i °i
-o s
I
1
P.
os
tu
o .1
.1.5
a = 5"
S. Kuroda/J. Wind Eng. Ind. Aerodyn. 67&68 (1997) 239-252 249
5. Static force coefficients
In the present computations, the drag coefficient is non-dimensionalized with
respect to H, the lift coefficient with respect to B and the pitching moment coefficient
with respect to B2. Here B is the girder width and H the girder depth (see Fig. 1). These
aerodynamic coefficients have not been corrected for the side-railings and the central
crash barrier. Fig. 7 shows the time series of computed drag, lift and pitching
moments. Fig. 8 shows the comparison of static force coefficients between the compu-
tation and the wind tunnel test 2. Computational static force coefficients were
calculated by time-averaging the instantaneous values over the interval of t = 20-40
(for z~= 10, t = 4(~80).
As a whole, the computed drag and pitching moment agreed well with the measure-
ment from the wind-tunnel experiment. Especially, as for the pitching moment, the
agreement between both is fairly well. The computed lifts are also in good agreement
with the experiment in the range of positive angles of attack (wind from below).
However, in the range of negative angles of attack (wind from above), the computed
lift coefficients were somewhat deflected from the experimental data.
The major cause of this discrepancy in the ~ negative cases can be considered as
follows. In the ~ negative cases, the computed flow was attached to the upper surface
of the girder as seen in Figs. 4 and 5. However, in the wind tunnel test, the flow cannot
co
i-
z 1.5
o
i1
u_ 1
i..iJ
O
o
uJ 0.5
on-
O
u_ 0
o
i
z -0.5
>–

©

i..u

i

0n-

I.tJ

-1

-1.5

a

o CD(cal)

~ CL(cal)

CM(cal)

5 -10 -5 0 5 10 15

ANGLE OF ATTACK (deg)

Fig. 8. Comparisonof computedstatic forcecoefficientswith experimentaldata 2.

250

S. Kuroda/J. Wind Eng. Ind. Aerodyn. 67&68 (1997) 239-252

F–

Z

Iii

t.l_

I.J_

U.l

©

o

L.I.I

on-

Ou_

o

Z

c~

0r r

LLI

0.,5

-0.5

-1.5 I

L

Ot _.=. 0 °

~llll,ll~,~, I,,,, i,i,iI,,,,Ibll,i

o I0 20 30 40 50 60 70 80

TIME

Fig. 9. Time histories of computed aerodynamic forces at 2 = 0′.

be supposed to have been attached to the upper surface, because the side railings and

crash barrier were included in the experimental model. The cause of the difference of

the lift coefficient between the computation and experiment in the cases of c~negative

may be the difference of the flow pattern along the upper surface between the two. On

the contrary, for ~ positive, the railing on the windward side must have been buried in

the vortices generated at the leading edge (see Fig. 4) in the wind tunnel test and had

no significant effects on the overall flow field. Consequently, static force coefficients

showed a good agreement between the computation (without railing) and the wind

tunnel test (with railing) in the e positive cases.

Finally, a frequency characteristic of the computed aerodynamic forces is referred

to cursorily. In Ref. 2, no description can be found of Strouhal number of the fixed

H4.1 section (Fig. 1). However, vertical vortex shedding response for the selected

section H9.1, which was almost identical to the H4.1 section, in smooth flow is found

in Ref. -2. The reduced velocity (V/(fB)) corresponding to onset of vortex induced

oscillations was around 0.9 -2.From this, the Strouhal number based on the girder

width for the H9.1 section is estimated to be around 1.1. For comparison, a Fourier

transform of the computed lift coefficient at ~ = 0° is presented. To perform a fre-

quency decomposition by using a long duration data, the flow simulation of ~ = 0°

was continued until t = 80. The transform was taken in the range of t = 2~80. Fig. 9

S. Kuroda/J. Wind Eng. Ind. Aerodyn. 67&68 (1997) 239-252 251

0.05~, ‘ ‘ i ….

0.O4

n”

I–

oI.U 0.03

13…

CO

r’r

LU

0

I ….

,

i i i

< 0.02 O.O1 0 012345 FREQUENCY ( f B/U ) Fig. 10. Power spectrum of computed lift coefficients at :~ = 0". shows the time histories of the computed aerodynamic forces and Fig. 10 the power spectrum distribution of fluctuating lift coefficients at ~ = 0°. Although there is a peak at a nondimensional frequency of 1.15, no particular predominant component can be found in the spectrum. It is probable that such a distribution which has no predomi- nant component is physical, because the cross section is shallow. Further investigation is needed to draw a conclusion. 6. Conclusions The flow over the box girder section model of the Great Belt East Bridge at a high Reynolds number was numerically simulated. Computed static force coefficients were compared with the measurements from the wind-tunnel experiment conducted at DMI. Computed drag and pitching moments agreed well with the experimental data. Computed lift also agreed well in the range of positive angles of attack, but is somewhat larger at negative angles of attack. This discrepancy in the blow-down case is thought to be caused mainly by the absence of equipment such as side railings and crash barriers in the computation. The predicted lift coefficients in the negative cases may be improved by considering the presence of the bridge equipments. 252 S. Kuroda/J. Wind Eng. lnd. Aerodyn. 67&68 (1997) 239-252 References 1 T. Miyata, H. Yamada, K. Kazama, On application of the direct flutter FEM analysis for long span bridges, in: Proc. 9th Int. Conf. on Wind Engineering, 1995, pp. 1030 1041. 2 T.A. Reinhold, M. Brinch, A. Damsgaard, Wind tunnel tests for the Great Belt link, in: Proc. Int. Syrup. on Aerodynamics of Large Bridge, 1992, pp. 255 267. 3 A. Larsen, Aerodynamic aspects of the final design of the 1624m suspension bridge across the Great Belt, J. Wind Eng. Ind. Aerodyn. 48 (1993) 261. 4 S.E. Rogers, D. Kwak, An upwind differencing scheme for the time accurate incompressible Navier Stokes Equations, AIAA J. 28 (2) (1990) 253. 5 M.M. Rai, Navier-Stokes simulations of blade-vortex interaction using high-order accurate upwind Schemes, A1AA Paper 87-0543, 1987. 6 J.L. Steger, D.S. Chaussee, Generation of Body-Fitted Coordinates using Hyperbolic Partial Differen- tial Equations, SIAM J. Sci. Stat. Comput. 1 (4) (1980).