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 –

©

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