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

Best services for writing your paper according to Trustpilot

Premium Partner
From $18.00 per page
4,8 / 5
Writers Experience
Recommended Service
From $13.90 per page
4,6 / 5
Writers Experience
From $20.00 per page
4,5 / 5
Writers Experience
* All Partners were chosen among 50+ writing services by our Customer Satisfaction Team

Shinichi Kuroda 1

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


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

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



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




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)


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


krlyp + t,V + rhVJ”


= – j*” (2b)




a=j ,


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)


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 >–







o CD(cal)
~ CL(cal)

5 -10 -5 0 5 10 15

Fig. 8. Comparisonof computedstatic forcecoefficientswith experimentaldata 2.


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






0r r



-1.5 I

Ot _.=. 0 °

~llll,ll~,~, I,,,, i,i,iI,,,,Ibll,i
o I0 20 30 40 50 60 70 80


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 ….


oI.U 0.03




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).