Computational fluid dynamics as an early aerodynamic design tool for a typical bluff bridge deck Open Access

$B$

section width

$CD$

drag force coefficient

$CL$

lift force coefficient

$CM$

moment coefficient

$c$

damping coefficient

$D$

section depth

$EF$

fluid energy

$Ek$

kinetic energy

$ENET$

net energy

$Ep$

potential energy

$ES$

system energy

$ETOT$

total energy

$E*$

reduced energy

$FD$

drag force

$FL$

lift force

$fc$

prescribed oscillation frequency

$fn$

natural frequency

$fvs$

vortex shedding frequency

$k$

stiffness

$L$

span-wise length

$MZ$

pitching moment

$m$

mass

$Q$

quality factor

$Re$

Reynolds number

$St$

Strouhal number

$t$

time

$tr$

reduced time

$U$

free-stream flow velocity

$Ur$

reduced velocity

$Y$

displacement amplitude

$Yr$

reduced displacement amplitude

$y$

displacement

$yr$

reduced displacement

$α$

angle of attack

$ρ$

air density

$ω$

natural angular frequency

For heavy and stiff bridges (e.g. Rialto Bridge in Venice, Italy or Tower Bridge in London, UK), static wind loads are of most concern, while for lighter and more flexible designs, dynamic wind loading effects become the limiting issue. Broadly, there are two mechanisms of concern: vortex shedding-induced vibration (VIV), which is caused by periodic aerodynamic forces resulting from wake instability; and coupled flow–structure instabilities, including galloping and flutter. In addition, slender bridges can be affected by turbulence-induced vibrations, that is buffeting.

While large-amplitude vibrations can be managed with retrofitted tuned mass dampers (TMDs) once a problem arises, a more common approach is to eliminate the aerodynamic susceptibility at the design stage. This is done by investigating the aerodynamic profile of the bridge as a whole, but also its elements individually (e.g. deck, pylon, cables). Since the infamous collapse of the Tacoma Narrows Bridge, structural engineers have been striving towards a better understanding of the highly non-linear bridge–wind interaction. Although many design codes and standards – for example, CD 363 (Highways England, 2020), EN 1991-1-4 (CEN, 2005) and PD 6688-1-4 (BSI, 2015) – give guidelines for the consideration of the aerodynamic effects on bridges, for long-span bridges, conducting wind tunnel tests is still a state-of-the-art approach to achieving aerodynamically acceptable designs. However, at the preliminary design stage, computational fluid dynamics (CFD) may provide valuable information on the bridge’s susceptibility to aerodynamic effects.

There are many examples of bridge designs for which aerodynamic effects were insufficiently considered. Depending on the type of effects that occurred, the outcomes varied from complete structural failure at the extreme, to temporary closures which were resolved with some form of retrofitting. Perhaps the most famous case, which instigated the development of bridge aerodynamics, is the Tacoma Narrows Bridge and its collapse in 1940, only 4 months after it was opened to the public. Many researchers tried to explain the reason for its collapse and the majority of these attempts were summarised by Larsen (2000). He gave a final verdict explaining that the vortices separating at the windward edge descended along the floor web and eventually synchronised with the torsional motion. This happened at a wind velocity of approximately 19 m/s. In the same paper, he states that even if this phenomenon had been avoided by altering the leading edge, the bridge would have experienced the two-degrees-of-freedom flutter at the velocity of 30 m/s, which is still not acceptable for today’s designs. In general, VIV will not cause a structural collapse of a bridge, however, it can cause fatigue and serviceability issues. In 1980, the Rio Niteroi Bridge in Brazil experienced VIV with excessive vertical deck motion. A comprehensive wind tunnel study determined that it was caused by VIV for wind velocities of approximately 14 m/s. The bridge remained closed until it was fitted with passive and active dynamic control devices (Battista and Pfeil, 2000). In 2010 a bridge over the River Volga, Russia, experienced high-amplitude vertical vibrations (approximately ±400 mm), which forced authorities to close the bridge. To suppress future vibrations, TMDs were designed and installed (Benicke and Butz, 2015). Other cases of VIV can be found in the literature: Tokyo Bay Bridge in Japan (Fujino and Yoshida, 2002) and The Great Belt East Bridge in Denmark (Schewe and Larsen, 1998) both experienced VIV.

As has been seen, proper consideration must be given to aerodynamics during the bridge design process to avoid unforeseen closures and retrofitting at the service stage. Modern wind tunnel tests can identify possible issues caused by the bridge’s high aerodynamic susceptibility but are costly and time consuming. Having prior knowledge of what to expect from these tests, before committing to them, can yield more accurate and faster testing campaigns. Design codes provide guidelines (see, e.g. CD 363 (Highways England, 2020)) for simplified desk-top calculations which produce critical wind velocities for different aerodynamic phenomena (VIVs, galloping etc.). These calculations always have a certain level of embedded empiricism as they combine a series of wind tunnel testing and some form of curve fitting to generalise the results for wider applicability. The generalisation of calculations leads to poor confidence for decks that have geometry significantly different from the geometry of the decks used for the development of the codes, necessitating more comprehensive wind tunnel testing.

Another approach for obtaining a preliminary insight into the behaviour of the bridge under wind conditions is by undertaking CFD analysis. CFD employs numerical methods for solving the non-linear second-order partial differential equations describing a fluid flow. In practice, in bridge engineering, due to the size of the structures and wind velocities occurring, the flows are almost always turbulent. When the flow is in the turbulent regime, the velocity and other flow variables vary in time chaotically and randomly. For an easier description of the flow field, the velocity field variable can be decomposed into a steady mean value $u¯$ and a fluctuating component $u′(t)$ resulting in $ut=u¯+u′(t)$⁠, termed Reynolds decomposition (Versteeg and Malalasekera, 2007). Once the Reynolds decomposition is applied to all flow variables, a resulting set of equations is called Reynolds-averaged Navier–Stokes (RANS) equations. There is a great number of books and papers that describe the derivation of the RANS equations, here it is only noted that when the flow variables are decomposed, the system of equations describing the fluid flow becomes underdetermined.

To be able to compute the turbulent flows using RANS equations, it is necessary to employ turbulence models which are used to close the underdetermined system of equations. The most used turbulence models are the k–ε model, the k–ω model and the shear stress transport (SST) model.

A more complex way of dealing with the turbulence is by numerically resolving the flow, that is solving the equations describing the flow. This approach in CFD is called direct numerical simulation (DNS). To be able to resolve the small turbulence scales, the grid size requirements increase dramatically when compared to simulations employing the turbulence models and depending on the flow fluctuations frequency range, the time step size is decreased. With this, the processing power needed for running these highly computationally demanding simulations renders most of today’s workstations unfit for the task.

To avoid resolving the smallest eddies, large eddy simulation (LES) uses a spatial filter to separate small from large eddies and numerically resolve only the later ones. Alternative simplifications are offered by delayed detached eddy simulation (DDES). DDES is a hybrid RANS–LES model which employs standard RANS turbulence modelling in the near-wall region while in other regions it enables LES. Although this approach decreases the computational power needed for the LES simulations, the domain should be modelled as three-dimensional (3D), instead of the two-dimensional (2D) approach usually employed with turbulence models. Two-dimensional LES simulations produce significantly worse results than their 3D counterparts when vortex shedding is studied (Murakami and Mochida, 1995).

Although CFD has advanced dramatically, both through the development of sophisticated solvers and turbulence models and through the development of hardware capable of solving non-linear partial differential equations, the uncertainties and computational requirements that come with it still put it out of reach for many non-research endeavours. Nevertheless, many studies have tried using CFD to replicate the results obtained in the wind tunnel tests of bridges. Zhang et al. (2021) modelled the Rose Fitzgerald Kennedy Bridge in Ireland using OpenFOAM. Their analysis consisted of developing 3D models and using the DDES method with the SST for turbulence modelling. The mesh they developed for the analysis comprised approximately 33 million cells and for running the analyses they used a supercomputer based in the Irish Centre for High-End Computing (ICHEC). They compared the k–ε, the Spalart–Allmaras (SAS) and the standard SST turbulence models and found that the results most similar to those from the wind tunnel tests were produced when employing the SST turbulence model. Furthermore, they used DES and compared its time-averaged results to the results produced using the SST turbulence model and found no large increase in the quality of the results, however, with the increase in computational time by a factor of 60. Nieto et al. (2009) studied the aerodynamic behaviour of a bridge deck at the conceptual design stage and provided the results for the static aerodynamic coefficients for a range of angles of attack. Based on that, they made recommendations on the most favourable shape to be used at the preliminary design stage (see Figure 1). The initial deck shape produced a lift coefficient with a negative slope around 0° angle of attack, which points to issues with galloping. They concluded that by adding baffle plates (vertical plates at the deck soffit) and streamlining the leading and trailing edges by adding fascia features, the final deck shape would produce a lift coefficient with a positive slope for the whole range of angles of attack and hence would be a better suited deck shape (from the aerodynamics perspective). These CFD models consisted of approximately 500 000 cells.

Zhang et al. (2020) used CFD analysis during the development of the VIV mitigation measures for the Yanpingba Yangtze River Bridge in China. Their models were 2D unsteady Reynolds-averaged Navier–Stokes (URANS) using the SST turbulence model. They used the CFD results to identify the regions of the deck responsible for the generation of the vortices and to better understand the structure of the vortices forming in the wake. This allowed them to develop several VIV mitigation measures, which were later implemented in the wind tunnel testing. All their simulations were done on a stationary deck model. Bruno and Mancini (2002) analysed the importance of barriers on the aerodynamics of two semi-streamlined bridge decks (see Figure 2) by using URANS (with the SST model) and LES. They showed that by adding the barrier, the drag coefficient had increased but also the mean lift coefficient was lowered; the sections prone to vortex shedding experienced a reduction in the amplitude of the lift fluctuations. They concluded that the numerical simulations could be a useful complementary method of aerodynamic analysis of bridge deck shapes.

Two-dimensional LES analysis was used by Bruno and Khris (2003) in their analysis of the Great Belt Bridge (Figure 3). They managed to obtain results comparable to those from the 3D LES analysis, however, they state that the methodology should not be applied to more complex flows before similar validation with the 3D LES is undertaken and that it might not yield results of the same quality.

Twin box deck sections were also subject to several CFD analyses. Goering and Ramponi (2019) used the SST turbulence model to identify flutter derivatives for galloping and static force coefficients for a range of angles of attack for the Stonecutters Bridge in Honk Kong (Figure 4). They conclude that ‘low-cost’ CFD models can provide design teams with significant information when making decisions on the bridge deck geometry. Jeong et al. (2019) used 2D URANS and the SST turbulence model to analyse a twin deck section (Figure 5) for angles of attack ranging between −10° and +10°. Their model comprised approximately 300 000 cells, and the stationary deck section. The results included static force coefficients and the Strouhal number estimate. It was concluded that, using the SST turbulence model, the CFD simulation was able to recreate most of the wind tunnel results with acceptable accuracy. Some discrepancy was observed at higher angles of attack, and this was attributed to the two-dimensionality of the CFD model.

Zhan et al. (2014) used the 2D LES model to analyse the fluid–structure interaction (FSI) problem of a twin box bridge deck (geometry similar to that of the Stonecutters Bridge in Figure 4). They obtained the results for the lock-in range and the displacement amplitude due to the VIV. The results were of the same order of magnitude as the ones from the wind tunnel tests (see Figure 6). The discrepancy between the CFD and the wind tunnel results was attributed to the existence of the non-linear damping in the wind tunnel model.

A similar analysis, again for a streamlined deck section, albeit with the barriers, was undertaken by Chen et al. (2021). They employed the SST model on a 2D grid comprising approximately 240 000 cells to estimate deck VIV response at various angles of attack. The models were used to inform on the barrier shape and size. It was found that the barrier plays an important role in the VIV response of the bridge deck. Increased porosity of the barrier increases the negative pressure region behind it. This means that the vortex that is the primary cause of vertical motion will transfer more energy to the structure. As a result, the bridge deck will be more susceptible to VIV issues.

It can be concluded that CFD has been used in bridge engineering in past decades; however, most of the analyses have been done on streamlined and semi-streamlined stationary deck sections (as seen in Figures 1–5). At a low angle of attack, the flow field around streamlined sections is less disturbed than in the case of bluff sections, which makes mesh generation an easier task. An example of the bluff deck section (Figure 1) is that analysed by Nieto et al. (2009), who reported that their model had approximately half a million cells, while the more streamlined sections (see e.g. Figure 5) were analysed by models comprising approximately 300 000 cells. Also, due to the complexity that the barriers introduce, most of the analyses omit them, while it was shown that the barriers play an important role in the aerodynamic behaviour of the bridge deck (Bruno and Mancini, 2002). Even fewer studies have been published on the coupled FSI problem of the bridge decks. While the static force coefficients provide useful information on the stability of the deck, the lock-in range and the amplitude of the displacements are what govern the bridge design when it comes to the VIV phenomenon. Furthermore, to the best of the authors’ knowledge, only Nieto et al. (2009) and Zhang et al. (2020) published their work on using CFD as a tool to inform the design at the early, conceptual stage.

Experience of practice suggests that the reasons for the lack of CFD presence at the design stage are twofold: (a) a detailed CFD analysis requires a large amount of computing power for it to be done in a time frame acceptable for the design stage; (b) the results achieved by CFD often bring uncertainties that can be removed only by committing to a wind tunnel testing campaign. Thus, engineers are still reluctant to choose CFD as a design tool in bridge engineering.

The aim of the present research is to develop recommendations for using CFD as a preliminary design tool for a typical bluff bridge deck in an industrial context. The models used here range from steady-state to dynamic FSI and prescribed motion models, each of them bringing an addition to the set of results that can be used to inform the design at an early stage.

The models must be of such size and complexity that they can be run on desktop workstations while requiring a reasonable amount of computational time, a constraint normally set by the industry rather than by academia.

Wind tunnel data used in the design process of the Northern Spire Bridge, over the River Wear in the UK, were available to identify and validate the most useful CFD modelling approaches. The bridge was selected because it comprises a deck that can be characterised as a bluff section. The wind tunnel campaign was executed in 2018 at Centre Scientifique et Technique du Bâtiment, Nantes, France in the boundary layer wind tunnel facility. At the time of the testing, it was not envisaged that the data would be utilised for research purposes, so the output of the tests is limited to what is required by the design and certification procedures (i.e. static force coefficients for a range of angles of attack and vibration amplitudes at a range of wind velocities). The data were presented in the factual report (Flamand and Knapp, 2018), which is not publicly available. Details are given here so that the experimental data can be better understood but these are not the final conclusions and recommendations of the factual report.

The tests were done at a scale of 1:50 at a Reynolds number of approximately 2 × 10⁴ based on the deck depth (⁠$D$⁠) with a velocity scale of 4.35 (defined as a product of frequency and geometric scales). The deck shape and the adopted sign convention are shown in Figure 7. For the sprung model, the system damping ratio was measured at 0.45% of the critical damping. The span-wise length of the model (⁠$L$⁠) was 1935 mm.

The tests were conducted in two stages: tests of the fixed model, which yielded aerodynamic force coefficients and gradients; and the ‘free vibration’ tests with the model supported by springs (as shown in Figure 8), to obtain the response to periodic vortex shedding. The dynamic model had two degrees of freedom, heaving or bending (translation in the $Y$-direction) and rocking or torsion (rotation about the $Z$-axis). The computational models developed in this paper neglect the latter degree of freedom as the heaving mode dominates the VIV of this structure (based on the results presented in the factual report).

The aerodynamic coefficients for a range of angles of attack at a specific wind speed are shown in Figure 9. Turbulence intensity at the inlet was measured at 4%. These results are broadly consistent with other similar studies found in the literature. For comparison, Goering and Ramponi (2019) and Zhu and Chen (2013) reported the lift coefficients in their studies to vary from −0.6 to 0.3 and −0.6 to 0.6; the drag coefficients ranged from 0.15 to 0.1 and 0.45 to 0.3, while the moment coefficients ranged from −0.1 to 0.1 in both studies.

Static force coefficients are defined in Equations 1–3. It should be noted that all three coefficients were normalised with respect to the bridge width (⁠$B$⁠).

where $CL$⁠, $CD$ and $CM$ are lift, drag and moment coefficients, respectively; $ρ$ is the air density; $U$ is the free-stream flow velocity; $FY$⁠, $FX$ and $MZ$ are forces in $Y$ and $X$ directions and the rocking moment around the $Z$-axis, respectively; $B$ is the deck width at model scale (0.488 m); $L$ is the model’s span-wise length (1.935 m).

It is straightforward to verify using the static coefficients in Figure 9 that the bridge deck will not experience galloping or coupled flutter in the velocity range under consideration. Therefore, the response of the sprung model to aerodynamic forces will be dominated by VIVs.

The level of turbulence intensity measured in the wind tunnel during this test was reported at 4%. The wind tunnel model dynamic properties, as reported in the factual report, are summarised in Table 1.

The experimentally obtained root mean square (RMS) of the amplitude response is shown in Figure 10. The values reported here are converted to full scale as is normal practice in bridge testing. The maximum displacement occurs at approximately 17 m/s, suggesting a lock-in range between approximately 14 m/s and 20 m/s. The Strouhal number (St) is estimated as 0.075. The VIV response of the wind tunnel model will be discussed in more detail below in Section 5.

It is worth reiterating that the purpose of the current study is to identify practical modelling choices to use CFD as a design tool for a typical bluff bridge deck. The approach taken is to compare the results obtained from a commercial RANS code with various modelling choices to the wind tunnel data described in Section 2. Throughout this study, three types of models were developed, as described below.

• ‘Fixed’ models were used to evaluate the modelling decisions made. Static force coefficients were compared to the results from the wind tunnel data described in Section 2.1. The fixed models were also used to estimate the Strouhal number for this deck.

• ‘Free-to-oscilate (denoted as FSI)’ models, which include FSI, were developed using the grid from the best-performing fixed model. A detailed description of these models is given in Section 3.5. They provided response amplitudes as the wind speed increased. In effect, these were numerical simulations of the wind tunnel tests conducted, so the direct comparison with the wind tunnel data is straightforward.

• ‘Prescribed motion (PM)’ models were used to further investigate the behaviour of the bridge deck and offer a faster way of estimating the VIV response at different damping levels and flow velocities.

Figure 11 gives an overview of all modelling choices made for the present study. Three choices for the geometry were considered: the bare deck (without the girders and barriers); the deck including the girders; and the full deck. The results for each are given in Section 4.1.1. Three options for the domain were explored, and the results are shown in Section 4.1.2. The ‘2.5D’ analyses represent a 3D simulation done on a narrow domain with periodic boundary conditions in the span-wise direction. It is noted that these simulations should not be mistaken for conventional 3D analyses. The purpose was to investigate whether large geometric features which deviate from the 2D deck profile (e.g. 3D features on the barrier, or cross-girders) would significantly modify the flow and disrupt the vortex shedding process. It should be noted that the turbulence model is still based on an isotropic assumption, and so spontaneous span-wise variation in the flow field would not arise without 3D geometric features. Three turbulence models are covered, with each bringing an extra level of robustness. The simplest one considered was k–ε, this model was then upgraded to SST, and the final turbulence model considered was transition SST (TSST), Fluent’s proprietary model. More details are given in Section 3.4 and results are summarised in Section 4.1.3. Finally, three different analysis types were considered: steady state (steady fluid flow), transient state (rigid structure but unsteady fluid flow) and full FSI (vibrating structure, unsteady fluid flow).

The idea was to develop a CFD model sufficiently robust to produce results comparable to the wind tunnel tests, but also simple enough so that it can be used in everyday industrial practice.

As seen in Figure 11, in total, 81 combinations can be made from the suggested choices, although not all were simulated. The selected models are listed in the simulation matrix (Table 2) and the results are presented and discussed in the following sections.

In the following sections, a detailed description of the grid used for model S03 is given. All other models were developed following a similar approach. The model dimensions adopted in this study correspond to the model dimensions used in the wind tunnel tests (see Figure 7).

While it is possible to use a single, fully unstructured, polyhedral grid, in the current study subdomains are defined and meshed with an unstructured mesh. There are two main reasons for this approach. First, it must facilitate a fully coupled FSI numerical model, albeit with a rigid deck on flexible supports. By defining a subdomain around the bridge deck, only that block would require remeshing or mesh deformation. Furthermore, the mesh block immediately adjacent to the bridge deck surface can be translated with the deck, so moving the remeshing to a less dense mesh block. Second, defining subdomains facilitates control of the localised mesh density, allowing mesh nodes to be concentrated in regions which are fluid-mechanically important (e.g. in the wake, or in the vicinity of the barriers). An example of the mesh obtained with this approach is shown in Figure 12.

A fully translating body mesh could be considered a better option since it allows for the mesh deformation to be avoided altogether. However, in the case when installation effects (i.e. ground topology) and atmospheric shear layer are to be considered, this approach would not be feasible.

Grid sensitivity was assessed to achieve a solution not dependent on the refinement level of the mesh. A summary of the results for model S03 is given in Table 3. As can be seen, the change in all static coefficients becomes less than 3% for the medium mesh refinement level and as such is adopted for all simulations discussed below. It is interesting to note that the moment coefficient appears to be more sensitive to mesh density. This is probably due to the barriers and the fact that the bridge deck width (⁠$B$⁠) is large compared to the depth (⁠$D$⁠), which is likely to be true for most bridge decks. The maximum $y+$ value for all simulations was below 1.0.

Dirichlet boundary conditions were used at all four sides of the domain with constant velocity set at the upper, lower and upstream sides and constant pressure at the downstream side. The flow velocity for all steady simulations was set at 5 m/s, resulting in a Reynolds number of the order of 2 × 10⁴. The angle of attack was controlled by setting two velocity components, in the global $X$ and $Y$ directions, as defined in Equations 4 and 5, respectively. The flow was uniform. All solid boundaries were modelled with no-slip conditions.

where $α$ is the angle of attack, $Ubx$ and $Uby$ are two velocity components, and $U$ is the free-stream velocity. The reason to consider non-zero angles of attack in bridge aerodynamics is two-fold. First, depending on the wider topology of the landscape, the bridge deck may experience wind shear and hence a downwash. Second, and more importantly, the galloping or flutter stability is determined by the gradients of the aerodynamic coefficients, widely known as the Den Hartog criterion (Den Hartog, 1957).

For most of the simulations, the transition SST (TSST) turbulence model was selected with smooth flow conditions at the inlet. The TSST is a proprietary model developed for Fluent. It combines four transport equations present in the classic SST model with an additional two transport equations: for intermittency and transition Reynolds number. More information on this turbulence model can be found in the Ansys Fluent user guide (Ansys, 2021).

Additional analyses were done using SST and k–ε turbulence models, and increased turbulent flow conditions in the upstream region.

As previously noted, the aim of the study is to develop a model that can be used employing modest computational power (i.e. a desktop workstation), hence the computationally expensive LES model was not explored.

The FSI was modelled by integrating the single-degree-of-freedom dynamic model (Figure 13) into Ansys Fluent. The model is defined by the equation of motion, Equation 6.

where $y$ is the vertical displacement of the body, and one and two overdots represent the velocity and acceleration of the body, respectively. The system mass, damping coefficient and stiffness are defined with $m$,$c$ and $k$⁠, respectively. The excitation force is denoted as $FL$⁠, which in the FSI model is the lift force.

As the built-in, six-degrees-of-freedom, free body motion feature in Fluent does not allow damping to be directly controlled, the motion was accounted for with an additional lumped parameter model. The lift force was obtained by integrating the pressure over the surfaces of the deck and two barriers. The equation of motion was solved for the vertical velocity of the deck using the Runge–Kutta fourth-order method. The rigid motion was then applied on the deck and barriers, and sub-grid levels L5 and L4 (see Figure 12). Remeshing was done only on sub-grid L3 during the motion of the system, while sub-grid levels L2 and L1 were kept static and rigid.

All dynamic properties reported in Table 1 were transformed to ‘per metre’ of the deck’s span-wise length and as such defined within the user-defined function (UDF). The UDFs were developed in such a way as to support parallel multi-core simulations.

The model developed for the FSI simulations was adapted to allow for the prescribed motion of the bridge deck. The motion is defined by applying a cosine function defined with Equation 7 to the deck velocity in the vertical direction at each time-step.

where $Y$ is the prescribed displacement amplitude (varied over a range of values); $f$ is the prescribed oscillation frequency (kept constant and equal to $fn$⁠, the first natural frequency of the model); and $t$ is time.

In this section, results from the simulations on the rigid models (S01–S42 in Table 2) are presented and discussed. Three sets of simulations were done in steady-state conditions and one in transient. The steady-state solution was obtained using the pseudo-transient method. It is considered that the steady-state solution has converged when the monitored aerodynamic coefficients have reached an asymptotic value, and the normalised residuals fall below the specified target (in all cases the threshold for the normalised residuals was set at 1 × 10⁻⁵).

First, the influence of the girder and the barriers on the static coefficients was investigated. It was shown how different numbers of holes in the barriers affect the aerodynamic profile of the deck while keeping the barrier porosity constant. Second, multi-dimension models were considered, and the three-dimensionality of the barrier was examined. The results from these, more computationally expensive models, were compared to the results from the 2D models. The third set of simulations included a study on different turbulence models available in Ansys Fluent. These included the k–ε model, SST model and TSST model. The final simulation on the static deck model was done under transient conditions, which provided an estimate of the Strouhal number.

A common strategy in any CFD analysis is to remove geometric features which are either complex or with small details, to reduce the overall mesh size and quality. To explore what level of geometric simplification is appropriate, the bridge deck as shown in Figure 7 is simplified by removing the girders and barriers, resulting in the section shown in Figure 14 (model S01) yielding an idealised 2D profile. Then the additional features (i.e. barriers above and girders below) are progressively added to determine whether they are necessary to achieve a reasonable comparison with experimental data.

The simulation done with this, much simplified, geometry shortens the meshing time and computational requirements dramatically. It in turn leads to shorter solution times and exhibits fewer convergence issues. Static aerodynamic coefficients are shown and compared to wind tunnel results in Figure 15(a). It should be noted that the wind tunnel results are for the full deck section (deck, girder and barriers). Even though there is a large difference in the geometry of the two models, a relatively good agreement between the wind tunnel lift and the simulated CFD lift is obtained, in terms of both the slope around the 0° angle of attack (which is important for stability) and the overall magnitude. As expected, the drag coefficient shows a large discrepancy in magnitude when compared to wind tunnel results. There are two reasons for this. First, the drag for both geometries is normalised with the bridge width (⁠$B$⁠), which is the same for the two cases. The total area projected in the $X$-direction for the bare deck is 71% smaller than in the case of the deck with the girders. This causes a factor of 3.4. It could be argued that the deck depth (⁠$D$⁠) should be used for obtaining the drag coefficient. However, in line with the common practice in bridge engineering, deck width is used for all three coefficients.

The second reason is that the barrier and girder will generate a proportionately larger wake. Despite the large discrepancy in drag and moment coefficients, this model can be useful for an initial estimation of the lift coefficient behaviour with varying angles of attack.

The girders are geometrically simple, so including them does not dramatically increase the mesh size (an increase of about 2% in the number of cells). Figure 15(b) shows static coefficients for the deck section without the barriers (model S02). As with the bare deck analysis, values of drag coefficient are lower than their experimental counterparts, although the discrepancy is much smaller. Considering that most meshing difficulties occur in the vicinity of the barriers, this can be regarded as another model worth exploring before moving to the full deck analysis.

Finally, a representation of the barriers was added with four openings (model S03). The results of this simulation are seen in Figure 15(c). The drag and moment coefficients now follow the trend of the experimental data well, although with a slightly lower magnitude. The lift, however, is not significantly better than that obtained for the bare deck, and in terms of both gradient and magnitude around $α=0°$⁠, it is arguably worse.

It is clear that the lift force is not very sensitive to the exclusion of either the barrier or the girders, while the experimental values are matching the trend exactly (i.e. gradient), and the magnitude is comparable. However, both the drag and the moment are progressively improved by the inclusion of first the girders and then the barrier. This demonstrates that the modelling of the bridge deck without the girders and barriers is not appropriate.

While this study only deals with a traffic/pedestrian barrier of relatively small height, it is worth noting that a much taller, wind-shielding barrier may be present on bridge decks, affecting the aerodynamic profile even more. Furthermore, wind fairings and guide vanes may also be found on bridge decks; these are commonly employed for improving the aerodynamic performance of bridge decks and, as such, their effects are not studied here.

As can be seen, increased geometric complexity and fidelity to the overall bridge deck significantly improves the accuracy of the simulated data. Including a representation of the barriers is necessary to achieve a reasonable estimate for drag and lift. However, the barriers in the wind tunnel model have a 3D pattern, as can be seen in Figure 16, which is itself an approximation of the full-scale barriers. The wind tunnel barriers are chosen to have an equivalent blockage ratio to the full-scale design. In the 2D approximation, the barrier is modelled as a solid fence with a number of openings.

The size of the openings is varied to assess whether the number of openings in the barrier makes a difference for any of the static coefficients when compared to the no-barrier case. The findings are presented in Figure 17 where the change in error with the introduction of the barriers is compared for six different configurations. As expected, the drag coefficient is least susceptible to the varying number of barrier openings since the total area over which the drag force is integrated is kept constant. It is also noted that the introduction of the barriers positively affected the drag coefficient estimate by decreasing the error for all angles of attack. This too suggests that by changing the number of openings, it is only possible to introduce a limited benefit in the accuracy of the drag coefficient.

Varying the barrier configuration has a dramatic effect on the flow adjacent to the top of the bridge deck due to jetting through the lower opening. Thus, the lift and the moment coefficients exhibit significant variations for the five barrier configurations considered. The models with fewer barrier openings produced better lift coefficients when compared to the wind tunnel data, while the opposite holds for the moment coefficients. Thus, there is no clear ‘best option’ when it comes to the number of openings in the barrier. The four openings model has been chosen as a compromise configuration as it also has the benefit of being geometrically similar to the wind tunnel model, which has between three and four openings over the depth (see Figure 16).

As already noted, the barrier has a 3D geometry. Furthermore, on the underside of the bridge deck, there are cross-girders which are stream-wise. Although the treatment of turbulence is still based on an isotropic closure, the small-scale features of the barrier and the secondary flows that would be expected in the junction of the cross-girders with the deck may affect the far wake and hence the overall aerodynamic behaviour. To assess the impact of the 2D assumption, two 3D domains were considered. The first, which is referred to as the ‘2.5D’ model (see Figure 18), represents a narrow strip of the deck with periodic boundary conditions on the two stream-wise domain boundaries across the deck width (model S21). The geometry is still substantially of a 2D character. The strip width is such that it captures the periodic shape of the barrier, but the cross-girders are ignored.

The second model is referred to as the ‘3D model’ (see Figure 19) and it is represented by a strip of such width that the half distance between two cross-girders is covered (model S22). Two planes perpendicular to the deck are set as symmetry boundary conditions. With this approach, the domain for the simulation was four times wider than the initially meshed one, see Figure 20(a). As seen in the figures, the barrier for the 3D models closely mimics the wind tunnel model.

For both models, a polyhedral mesh is used as it decreases the number of cells required to achieve an acceptable solution, see Figure 20(b). The mesh sensitivity analysis for this mesh has not been done due to high computational requirements. However, the cell size chosen for these models is comparable to the cell sizes employed in the 2D domains.

Simulations presented here employed the SST turbulence model only (as shown in the simulation matrix in Table 2).

The results of the ‘2.5D’ model (model S21) and the 3D model (model S22) are shown in Figures 21 and 22, respectively. Both are compared to the wind tunnel data and the results from the 2D analysis (model S32).

There is not a substantial difference between the multi-dimensional models (‘2.5D’ and 3D) and the 2D models, suggesting that the extra computational expense is not justified; each simulation took approximately 30 CPU hours using 120 CPU cores. However, it should be noted that turbulence is treated with an isotropic model, so the highly 3D turbulent wake will be poorly captured. A full DDES or LES model may yield improved results (Zhang et al., 2021), but the computational effort will increase dramatically, reducing the number of configurations that can be investigated. It is also noted that should a 3D model be of larger span-wise length than the ones covered here, the findings could be different. However, the purpose of this task was to investigate whether a limited extrusion to the third dimension could yield results better than those obtained from the 2D analyses.

To understand why the increased geometric complexity does not yield improved aerodynamic coefficients, the flow field of the near wake has been visualised. Two shear layers, coming from the top and the bottom of the deck, encapsulate the flow features originating from the barrier and limit their effect to the region just downstream of the deck, while the far wake is unaffected. This is substantiated in Figure 23 where the $Y$-vorticity contour for the plane at 2.6 m from the deck centroid and the $Y$-velocity streamlines are shown.

The simulation approach taken is steady RANS, as this is substantially less computationally expensive than LES or DDES. However, it is well known that the choice of the turbulence model in RANS may have a significant impact on the results. The models considered here are:

• the k–ε model S31

• SST — model S32

• TSST – model S03

TSST is a proprietary model developed by Fluent, see Ansys (2021) for more details. It should be noted that the Reynolds number of the wind tunnel (based on deck depth, $D$⁠) is in the transitional range at the lower velocities (∼10⁴), but the Reynolds number regime for the full-scale bridge will be fully turbulent (∼10⁶). Low Reynolds wall treatment is employed in all simulations. Figure 24 shows static coefficients for a range of angles of attack for different turbulence models employed and compared to the wind tunnel data. As can be seen, the TSST model yields the closest result compared to the wind tunnel data, especially for the lift coefficient.

While the SST model is developed for fully turbulent flows (⁠$Re>105$⁠), the TSST model outperforms it when the flow is in the transition region, between the laminar and fully turbulent conditions. This is documented by Lodefier et al. (2004).

Bridges operate in the atmospheric boundary layer, and so the upstream airflow typically exhibits high levels of turbulence. In the wind tunnel, this is modelled with turbulence generators. Figure 25 shows the static coefficients obtained with the TSST turbulence model for a turbulence intensity (⁠$Iu$⁠) of 10% at the inlet, which resulted in the $Iu$ of approximately 4% just upwind from the deck (model S41). The dissipation of turbulent kinetic energy in the turbulence model tends to reduce the free-stream turbulence levels more rapidly than is experienced in practice. In the same figure, static coefficients for the smooth flow (⁠$Iu=0%$⁠) at the inlet are shown (model S03), and a clear difference between the two models is noticeable. Both the magnitude and the slope of the curves show better correspondence with the wind tunnel data for a turbulent inflow.

Similar observations can be made for the simulations of the 3D SST model with increased turbulence intensity (model 52) shown in Figure 26. Here the change in lift slope around 0° angle of attack is even more noticeable and important. Having the slope of the lift coefficient negative around 0° angle of attack would lead one to believe that the test may be susceptible to instability.

The mesh and geometry of the fixed model (model S03), which was used to obtain the aerodynamic coefficients using steady flow, can now be used to easily yield basic information about the vortex shedding by conducting transient simulations. The initial time-step size was gradually decreased until a constant value of the Strouhal number was reached. Figure 27 shows the time history of the static force coefficients obtained from model T01 where the time is reduced with the free-stream flow velocity and the deck depth (⁠$D$⁠). The frequency analysis of the lift coefficient signal is shown in Figure 28. It should be noted that only data after the reduced time of 12 are used (Figure 27). This is to avoid artefacts associated with the transient part of the dataset associated with the starting flow. At this time, the wake extends approximately ten deck lengths downstream. The Strouhal number (⁠$St$⁠) is defined as in Equation 8, where $fvs$ is the shedding frequency, $U$ is the free-stream flow velocity and $D$ is the deck depth. The transient analysis of the fixed deck model estimates $St$ at 0.065, which compares well with that obtained from the wind tunnel data shown in Section 2.

The mesh and geometry identified with the fixed model are now used as the basis for a FSI model, where the bridge deck is supported on a flexible structure as described in Section 3.4 above.

The deck depth $D$ and the system’s natural frequency $fn$ were used for the reduction of all variables within this section. Reynolds number was calculated with respect to $D$⁠. Reduced variables are defined in Equations 9, where $Yr$⁠, $Ur$ and $tr$ are reduced deck displacement, reduced flow velocity and reduced time, respectively.

Two turbulence models were used for the FSI models: SST (model D01) and TSST (model D02). Figure 29 shows the deck response obtained from these simulations and compared to the wind tunnel data. The SST model provided good agreement in the shape of the lock-in region. The TSST turbulence model performed considerably better than the SST in determining the velocity at which the lock-in occurs. The vibrations start to be noticeable at a reduced velocity of approximately 13.5, which is in good agreement with what is reported in the wind tunnel data. However, with increased flow velocity and for the value of Reynold’s number of approximately 13 × 10³, the oscillation amplitude starts to diverge from the wind tunnel values. Clearly, at these large amplitudes, the simulation is either overestimating the excitation or underestimating the dissipation in the flow. This may be due to the 2D assumption compounded by the simple isotropic turbulence model, as at large amplitudes the span-wise correlation will become more important.

The lock-in region for these two simulations is further examined in Figure 30. It shows reduced shedding frequency ( $fvs/fn$⁠) and the Strouhal law for this deck section. The Strouhal law, as defined in Equation 8, is based on $St=0.065$⁠, as obtained from the transient model (model T01). A noticeable difference is observed between the two models. The SST model predicts the shedding frequency locks to the system’s natural frequency before the coincidence, while TSST shows behaviour more similar to what is observed in the wind tunnel with lock-in starting at coincidence between the Strouhal and the natural frequencies. Figure 31 shows the spectrogram of the deck displacement obtained in the wind tunnel. For consistency, all the values are reduced as previously described and the spectrogram in Figure 31 is shown for the full range of velocities tested in the experiment. It can be seen that the lock-in occurs at the reduced velocity of approximately 13.5, which is more or less the coincidence between the Strouhal frequency and the natural frequency. The lock-in range, as indicated by the suppression of the Strouhal frequency, extends to a reduced velocity of approximately 20. Note that immediately before coincidence, the vortex shedding frequency is distinct from the natural frequency. The data from the TSST (model D02) in Figure 30 shows similar behaviour in both pre- and post-coincidence. The SST results, in contrast, show strong pre-coincidence lock-in, which is not supported by the wind tunnel data.

In Figure 31, at a reduced velocity of 7.5, there appears to be another small lock-in zone. As the Strouhal frequency at this velocity is exactly half the natural frequency, it is suspected that this is due to the fluctuating drag, which will be at twice the shedding frequency and hence will coincide with the natural frequency.

The difference in behaviour for the two flow regimes (locked-in and locked-out) is also illustrated in Figure 32, depicting the time record of the lift coefficient and the reduced deck displacement. Figure 33 shows the frequency domains for these two signals. Figure 32(a) shows a much cleaner response with two signals synchronised. This is confirmed in the frequency domain analysis shown in Figure 33(a). It should be noted the purpose of the power spectral density (PSD) is to show the difference in peak frequency between displacement and vortex shedding excitation. The effect of spatial and temporal quantisation can be seen mainly in the higher frequencies, while the effect of non-stationarity caused by the relatively short time window will be seen primarily in the lower frequencies (this truncation error will also cause a reduction in the peak values of the PSD). Thus, while the absolute amplitude of the PSD will be changed, the apparent location of the dominant frequency will not be affected.

Figures 34 and 35 show the vorticity in the flow field at lock-in (reduced velocity of 17) and lock-out (reduced velocity of 20), respectively. In both cases vortex shedding in the wake is evident, but in the lock-in case, it shows spatial periodicity as would be expected. When the flow is in the locked-in regime (Figure 34), large alternating vortices form in the near-deck wake, within a deck width (⁠$B$⁠). The vortex street, while very ordered, exhibits what looks like a 2P mode on the lower side, but a 2S mode on the upper vortices. In the case of locked-out flow (Figure 35), the vortices are formed much later downstream and the structures never form a clear vortex street as the two frequencies (vortex shedding and structural natural frequency) are distinct and incommensurate.

Although free-to-oscillate simulations provided valuable information on the maximum displacement amplitude and lock-in region under VIV, the simulations must be performed in series. That is, due to the nature of the lock-in phenomenon, the analyses must be initiated at a velocity lower than the expected onset velocity. Because of this, the time required to reach the maximum displacement amplitude can be considerable. Many studies have employed prescribed motion simulations to study the VIV of structures submerged in the moving fluid (see e.g. Jiao and Wu, 2018; Paidoussis et al., 2014; Staubli, 1983). The two main advantages of these simulations are that they can be run simultaneously and that the important properties, such as displacement amplitude and oscillation frequency, are defined (i.e. can be controlled).

When the limit cycle oscillations (LCOs) are reached, the net energy of the dynamic system during one cycle equals zero. That is, during one oscillation cycle, the energy dissipated by the structure and the work done on the body by the fluid forces are equal.

In this section, the prescribed motion model (model D03) is used to evaluate the work done on the bridge deck by the vortex-induced forces. This is then compared to the dissipated energy that originates, primarily but not only, from the structural damping of the dynamic system.

The total energy of a damped harmonic oscillator is defined by a sum of the potential and the kinetic energy. Equations 10–12 define potential, kinetic and total energy, respectively.

where $γ=c/m$ is the damping–mass ratio and $ωd$ is the natural frequency of the damped system. For lightly damped systems, $γ<<ω0$⁠, the cosine in Equation 12 can be neglected. This reduces the total energy to Equation 13.

Introducing the quality factor (Q factor) as $Q=ωd/γ≈ω0/γ$ and evaluating the energy change during one oscillation period, $T=2πω0$⁠, the total energy loss during one oscillation cycle is given by Equation 14.

The energy brought into the system by the work done on the bridge deck by the vortex shedding forces during one oscillation cycle can be evaluated by integrating the lift power over one oscillation cycle, as given in Equation 15.

where $y˙$ and $FL$ are the deck velocity and fluid forces exerted on the deck in the vertical direction, respectively.

Finally, the total net energy in one oscillation cycle is given in Equation 16; it comprises the energy brought in by the fluid (⁠$ΔEF$⁠) and the energy taken away by the structural damper (⁠$ΔES$⁠).

For an easier comparison and discussion, the energy is non-dimensionalised as in Equation 17.

where $U$ is the wind velocity; $ρ$ is the air density; $D$ and $B$ are the depth and width of the bridge deck; and $Y$ and $f$ are oscillation amplitude and frequency, respectively.

In the current study, it is assumed that at lock-in (which is the region of primary interest), the ratio between the shedding frequency and deck oscillation frequency is unity. Many papers confirm that this is the case for high mass ratio systems (see e.g. Williamson and Govardhan, 2004). It was also assumed that the fluid-added mass is negligible when compared to the mass of the deck. Equations 10–14 also assume that the motion of the body is purely sinusoidal and that higher frequency oscillations are negligible. This is a common assumption and is supported by Figure 28, which only exhibits a single frequency.

Figure 36 shows the aerodynamic energy transferred into the system every oscillation cycle, Equation 15, for a range of oscillation amplitudes and at the reduced velocity of 16.8. In the same figure, energy taken away by the damping, Equation 14, is shown with the dashed line. For a given damping level, marked in red is the unstable region, that is the region in which LCOs can occur. The green colour marks the region in which the energy taken out from the system is greater than the aerodynamic energy brought into the system, and thus the LCO cannot occur. In the same figure, aerodynamic energy calculated from the FSI system (model D02) for $Ur$ equal to 16.8 is shown. In theory, this point should coincide with the point where ‘positive’ and ‘negative’ energy contributions are equal. The discrepancy can be attributed to the assumptions embedded in Equation 14 and the prescribed motion model being defined to model only the first mode of oscillation.

A different way of presenting the same data is given in Figure 37. Here, the ‘positive’ and ‘negative’ energy contributions are subtracted defining the ‘stable’ and ‘unstable’ regions above and below the zero-energy line. Again, the FSI estimate of the energy should lie on the zero axis; however, because of the limitations previously discussed, this is not the case. An interesting observation can be made for $YrRMS<3.5%$⁠, the LCO does not exist for the analysed flow velocity and the given damping level. So, to enter the LCO, the system needs to have an initial displacement, large enough so that the system enters the ‘unstable’ region. With this, the aerodynamic energy per cycle is greater than the dissipated energy, and the oscillation amplitude increases until the equilibrium between aerodynamic and dissipated energy has been reached. This happens at approximately $YrRMS=8.5%$⁠, which is the displacement amplitude estimated by the FSI model.

Here it has been shown that, for a given flow velocity and damping level, it is possible to predict, relatively accurately, the maximum oscillation amplitude by using the prescribed motion model (when compared to the FSI model). The main advantage of this approach over the FSI model is that, for each data point, the simulations can be run independently and at the same time. This shortens the time needed to generate an estimate for the maximum amplitude. However, when it is unknown at what velocity the maximum amplitude occurs, one is required to run the simulations over a range of velocities and amplitudes generating a 3D representation of the aerodynamic energy that is transferred into the system.

Figure 38 shows the data simulated for a range of $Ur$ and $Yr$⁠. The aerodynamic energy brought into the system every oscillation cycle is shown in green, while the negative energy, taken out from the system during each oscillation cycle by the damping, is shown in red for a damping ratio of 0.45%. Where these two surfaces intersect, the system energy summed over one oscillation cycle equals zero and the LCO occurs. The same data, shown from the top view, are given in Figure 38(b). The region where the aerodynamic work is larger than the negative work done by the damper is marked in green. The opposite is true in the region marked in red. In addition, the results from the FSI model (model D02) are plotted in black. It is seen that the intersection between the red and green regions, where the ‘positive’ and ‘negative’ energy flux are equal and the LCO occurs, matches closely the FSI model predictions of the response amplitude. This approach also predicts the oscillation amplitudes which exist at the LCO for each flow velocity, and they are in good agreement with the FSI model.

The main advantage of the energy approach is that the effect of different structural damping can be easily and rapidly evaluated. For example, in Figure 38(a), the blue surface represents the energy dissipated with a damping ratio of 0.55%. The projection onto the velocity-response amplitude plane can be seen in Figure 38(c), which shows a reduced peak response amplitude, and a slightly reduced lock-in range, as would be expected for increased damping.

Figure 39 shows the aerodynamic energy brought into the system with every oscillation cycle for $YrRMS$equal to 7% and 11% at $Ur=16.8$⁠. For the first case, the energy value stabilises after approximately 25 cycles and remains positive, meaning the fluid force is exerting work on the deck. The opposite holds for the second case, with $YrRMS=11%$⁠, the energy brought into the system is negative, meaning that the deck is exerting work on the surrounding fluid. Having negative aerodynamic work excludes this ‘$Ur$ – $YrRMS$’ combination as a possible scenario under which LCO can occur, regardless of the structural damping. This is confirmed in Figure 36. To have the curve showing the work done by the system damper intersect the curve describing the aerodynamic work, the system damping would have to be negative. For a given flow velocity, it is possible to estimate the maximum LCO amplitude that theoretically exists when the system damping is zero.

The results from the prescribed motion analysis, for a given flow velocity, can be used to estimate what is the necessary system damping at which the displacement amplitude is acceptable. The flow velocity at which the maximum displacement occurs can be estimated based on the Strouhal number acquired from the static transient analysis. The transient analysis (model T01) gives an estimate of $St = 0.065$⁠, which translates to $Ur=15.4$at the coincidence. However, the largest displacement occurred at $Ur=16.8$⁠. Looking at the results from the FSI model (Figure 29), it can be noticed that after the coincidence, the shedding frequency remains locked to the natural frequency of the bridge, until it locks out at approximately $Ur=17$⁠. This means that, after coincidence, the increase in displacement comes from the increase in aerodynamic energy, which is increased due to a rise in the flow velocity. This behaviour continues until the shedding frequency locks out, which disrupts the inflow of the aerodynamic energy and the oscillations due to vortex shedding disappear. This sudden lock-out is noticeable in the FSI simulations and the prescribed motion models.

The overall goal of the current study is to show different ways in which CFD can be used as a preliminary design tool to provide early aerodynamic results for a typical bluff bridge deck. At the preliminary stage, the design is not concerned with the maximum displacement amplitude caused by the VIV, but with the overall aerodynamic stability of the bridge deck. If the value of the maximum acceptable displacement amplitude is set, the prescribed motion model can be used to evaluate if and for what wind velocity range that value will be exceeded. The maximum wind velocity up to which that range spans is set as the design wind velocity for the VIV defined by the design code. Figure 40 shows the aerodynamic energy per oscillation cycle for a set displacement amplitude over a range of flow velocities. The aerodynamic energy (⁠$ΔEF$⁠) for three reduced displacement amplitudes is shown with solid lines. They are indexed 1, 2 and 3, for $YrRMS$ values of 6%, 7% and 10%, respectively. The negative work done by the damping (⁠$ΔES$⁠) is shown with dashed lines for each of the three displacement amplitudes. The range of flow velocities for which the LCOs can occur decreases with the increasing displacement amplitude. For the highest displacement amplitude shown (⁠$Yr,3RMS=10%$⁠), the LCO cannot occur for the analysed range of reduced velocities, because the energy taken out of the system by the damper is always higher than the aerodynamic energy brought into the system by the fluid forces. This gives an indication of how the analysed deck will perform in the wind tunnel tests. As previously mentioned, the prescribed motion models can be run simultaneously for a number of data points so the graph shown in Figure 40 can be produced relatively quickly. However, it should be noted that in cases when the damping is rigidly set, and the design needs the maximum displacement amplitude due to VIV at an early stage, the FSI model will provide results faster.

In Section 1, it was shown that most of the CFD analyses done to date were undertaken on streamlined or semi-streamlined bridge decks. It has also been shown that only one provided input for the design process, while the others were part of wider research studies. The goal of this paper was to explore the applicability of CFD as a tool, which could be employed using a desktop workstation (i.e. modest computational power), in everyday bridge design to provide preliminary aerodynamic results for bridge decks whose shape is not aerodynamically optimised (i.e. bluff bridge decks). This will allow for more informed, more efficient and faster wind tunnel campaigns for these types of decks. A proposed methodology for using CFD in bridge design (at the preliminary stage) is given in Table 4 and summarised in Figure 41.

The key findings of this study are summarised below.

• Simplified deck geometry can produce useful initial results on the lift coefficient. However, the barriers and girders must be included in the analysis because the flow around the deck is heavily influenced by these elements.

• If the overall porosity of the modelled barriers follows that of the bridge design, the number of openings in the barrier can be varied with little effect on the static force coefficients.

• Expanding the 2D model to the third dimension only slightly changes the force coefficients while dramatically increasing the computational power needed. It was shown that two shear layers constrain the flow features originating from small 3D features in the downwind barrier, limiting the effect on the far wake. Thus, when turbulence is modelled and not resolved, 2D analysis can provide sufficiently accurate results. It is noted that, should a 3D model be of larger span-wise length than the ones covered in this study, the findings could be different, but LES or an anisotropic turbulence model may be necessary.

• For the simulations done at the Reynolds number of a similar order of magnitude to the one present in the wind tunnel test (∼10⁴), the TSST turbulence model provides the results closest to the wind tunnel tests. With increasing $Re$number (>10⁵), the TSST model converges towards the SST model. It should be noted that the TSST model requires increased computational effort due to two additional transport equations when compared to the SST. It was observed that, on average, the TSST models took approximately 15% longer to converge.

• A free-to-oscillate model (FSI) can provide results comparable to the ones obtained from the wind tunnel tests, successfully capturing the lock-in phenomenon and the range of velocities under which VIV occurs. The amplitude of the displacements is larger than what is seen in the wind tunnel tests. It is suspected that the overestimation of the excitation is due to the two-dimensionality of the simulations, which imposes a full span-wise correlation. Better results could be expected from a 3D analysis coupled with LES or an anisotropic turbulence model.

• The PM model, coupled with the energy balance approach (see Section 5.1), can provide results comparable to those obtained with the FSI model. The PM model produces results, per data point, faster than the FSI model. From the simulations done in this study, it was seen that, on average, the time needed to reach a steady solution for a single data point is halved. However, when the wind velocity at which the maximum displacement amplitude occurs is unknown, more data points are needed, making the FSI approach faster in yielding results. Nonetheless, the PM model is superior to the FSI model when the sensitivity to structural damping is to be investigated. It was also concluded that for a set displacement amplitude and damping level, the PM model can yield a range of flow velocities for which the VIV can occur.

To further test the applicability of the findings described above, the authors have recently obtained preliminary aerodynamic results for a separate bridge design with a deck of a similar shape to the one studied here. The analysis was done before the wind tunnel tests and the results were used for the development of the testing specification. As the design of this bridge is still in progress, the details of these analyses are not yet publicly available. However, the experience suggests that the approach is beneficial in the design cycle.

This paper demonstrates that CFD can be used for acquiring preliminary aerodynamic results in bridge engineering when used on bluff bridge decks of comparable shape to the one studied here. Further work should be undertaken to develop and investigate similar approaches that can be applied to a broader range of deck geometries.

The authors would like to thank Olivier Flamand of Centre Scientifique et Technique du Bâtiment for sharing the data from the wind tunnel analyses conducted on the Northern Spire Bridge model and promptly answering questions related to the wind testing itself.

The authors would like to acknowledge the Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support. The authors would also like to acknowledge Sunderland City Council and Farrans Victor Buyck CJV for their contribution to the paper.