1 Introduction

Long-span bridges are an essential component of modern transportation networks. For instance, the Fourth Road Bridge, which located next to the Queensferry Crossing, carried traffic on annual average daily of around 70,000 vehicles based on 2004/2008 inspections. The bridge was even closed for repairs and refurbishment to eliminate concerns over the lifespan [1]. The functioning of long-span bridges are, however, susceptible to adverse weather conditions, particularly high wind speeds which are unpredictable in nature. Wind is a critical source of bridge vibration as long-span bridges are slender and vulnerable to wind-induced vibrations [2]. A further complication is that the presence of a bridge (or any structure) can alter the surrounding wind environment. A recent example of wind conditions having an adverse effect on a long-span bridge is that of the Humen Pearl Bridge, China [3] with a maximum span of 888 m. The vertical amplitude of the bridge deck oscillations induced by wind induced vortices was found to be greater than 2.2 m which led to a complete shutdown of the bridge, and a requirement for further safety assessments. Due to the increasing impact of climate change such as uncertainties of regional wind speed variations and trends of storm intensity increasing [4], bridge aerodynamic studies have been attracting increasing attention recently in the design and operation of modern long-span bridges. In the meantime, the rapidly developed numerical weather prediction systems such as the Weather Research and Forecasting (WRF) Model are able to provide reliable wind conditions of cites, which can be a strong support for bridge aerodynamic studies.

Traditionally, wind tunnel tests are the preferred method of conducting bridge aerodynamic analysis. Scaled bridge models are tested in wind tunnels that have sensors installed to capture the wind flow data. This assists scientists and engineers in investigating aerodynamic features and wind field characteristics of bridges [5,6,7]. However, due to the significant cost of conducting wind tunnel tests and limitations such as effects from intrusive sensing equipment and scaling effects [8,9,10], scientists and engineers are keen to develop a more convenient, cost effective and accurate method of aerodynamic analysis.

The relatively recent and rapid development of computational capacity has supported researchers in exploring the possibility and effectiveness of performing numerical aerodynamic simulations rather than wind tunnel tests. Combining the fundamental principles of applied mathematics and fluid mechanics, computational fluid dynamics (CFD) has been recognized as an emerging method with great potential to model wind as fluid flow and to carry out detailed analysis of bridge aerodynamics. Two-dimensional (2D) CFD simulation studies of bridges can be traced back to the 1990s [11]. Various numerical studies [12,13,14,15] on circular cylinders have been able to capture aerodynamic phenomenon such as the drag force abruptly changing as the Reynolds number changes as well as vortex instabilities in the wake of a circular cylinder. The effects of wind fairing angles and bridge flutter derivatives were studied by comparing the results of CFD simulations of a scaled bridge deck section with wind tunnel test results [16,17,18]. A significant number of 2D simulations of bridges have been compared with wind tunnel experiments, and limitations of 2D simulations have been identified such as the neglection of effects of flow development and turbulence interaction in the spanwise direction. As a result of such limitations, researchers are now starting to consider the use of three-dimensional (3D) simulations to ascertain more accurately bridge aerodynamic behaviors [8, 19, 20].

To date, only a small number of 3D CFD simulation studies have been undertaken on long-span bridges. Most of them are scaled models that are designed to be compared with results of wind tunnel test reports and lack field data validations. Zhang et al. [21] conducted a 3D CFD model of the bridge deck of the Rose Fitzgerald Kennedy Bridge. Time-averaged aerodynamic coefficients of three different bridge decks were determined and showed general agreement with wind tunnel tests. However, some discrepancies existed which suggested that certain conditions of scaled model tests, such as blockage ratios at large angles of attack and sensor interference in wind tunnel tests, can cause inaccuracy of test results. Zhu et al. [22] created a 3D CFD model of the Queensferry Crossing Bridge deck, containing wind shields and sample vehicles. The wind effects from a range of yaw wind angles were considered. Subsequently, the aerodynamic coefficients of the vehicles were determined and compared with wind tunnel tests. It was found that elaborate mesh schemes with average y+ value lower than 4 can significantly support CFD simulations in accurately determine the aerodynamic coefficients. Tang et al. [23] modeled a 3D scaled twin-box bridge girder to analyze its flutter performance at large angles of attack. The authors found that, at a null angle of attack, the critical flutter wind speed increases and the flutter frequency decreases with the increase in width of central slot of the bridge deck. At large angles of attack, the upstream box tends to be the main source of torsional flutter instability as the central slot would divide large vortices into a smaller and more complex flow field. However, the analyzed flow field characteristics were performed in a scaled simulation model with results compared to wind tunnel test data. The authors noted that further validation with on-site monitoring data would strengthen their conclusions. Alvarez et al. [24] conducted 3D CFD simulations to gain some insights regarding the characteristics of the complex flow features around twin-box decks. The authors observed that the presence of the transversal connecting beam had a positive effect on decreasing the correlation of the force coefficient and spanwise length of the twin decks. Moreover, it was found that different slot cells can influence that correlation as well, which deserves further investigation. Kim and Yhim [25] performed a fluid structure interaction buffeting analysis on a 3D long-span cable-stayed bridge structure including the calculation of the dynamic responses and detailed wind flows. The properties of wind load distribution and abnormal aerodynamic phenomena can be determined. However, the method was computational expensive and it still required further verification.

The availability of field data for validation of models is a key shortcoming in all of the aforementioned 3D CFD bridge studies. In each case, the numerical models have been compared with results of scaled wind-tunnel test results. To the author’s knowledge, bridge aerodynamic studies incorporating full-scale 3D CFD simulations validated against field data have not previously been undertaken. In this paper, for the first time 3D full-scale CFD models for the middle span and central tower of the long-span Queensferry Crossing Bridge were developed and simulations are performed with the k − ω SST turbulence model. Raw field data collected from the Structural Health Monitoring System (SHMS) on the Queensferry Crossing was restored, selected and processed in the database management system PostgreSQL. The purpose of this work was to create and validate the 3D CFD model with field data to demonstrate that the full-scale 3D CFD model can accurately describe the true field flow conditions. WRF simulation results were used to provide the initial wind velocity conditions. A mesh sensitivity study and a cylindrical domain study were performed. The cylindrical domain shows good reliability and is much more convenient than common rectangle domain in studies that involves yaw angles. Simulation results were validated with the field monitoring data. The simulation results show good agreement with the SHMS field data which indicated that the model accurately captures 3D effects and the WRF simulation provides good estimates for wind velocity conditions.

2 Field data: the Queensferry Crossing

The long-span Queensferry Crossing shown in Fig. 1 (formerly the Forth Replacement Crossing) carries the M90 motorway which connects the traffic across the Firth of Forth between Edinburgh, at South Queensferry, and Fife, at North Queensferry, in the United Kingdom. The bridge was opened in August 2017 alongside the existing Forth Road Bridge. The 2.7 km total length makes it the longest three-tower cable-stayed bridge in the world. Each tower is 210 m making them the highest in the United Kingdom. The balanced cantilevers which extend 322 m north and south from the central tower are the world’s largest free-standing balanced cantilevers. The bridge deck carries two lanes of traffic in each direction alongside hard shoulders.

Fig. 1
figure 1

View of the Queensferry Crossing Bridge [26]

Full size image

Transport Scotland appointed Strainstall (later rebranded James Fisher Testing Services) to design and install a SHMS on the Queensferry Crossing to measure environmental actions, bridge actions and responses such as air temperature, relative humidity, strains, and displacements. In total, there are 1972 sensors of various types installed over the entire structure. There is an off-site monitoring system installed which enables the data to be managed remotely.

2.1 Anemometer sensors

There are 11 3D heated ultrasonic anemometers installed on the bridge structure (refer to Fig. 2). Three pairs of emitters and receivers are combined in a 3D anemometer to resolve the complete wind vector. These devices have a digital output that provides data on wind speed, wind direction and temperature per anemometer. The anemometers make use of ultrasonic sound waves to determine instantaneous wind speed by measuring how much sound waves traveling between a pair of transducers are sped up or slowed down by the effect of the wind. The speed of sound is dependent on the density of the air, which is principally determined by air temperature. Therefore, the sonic measurement can provide an estimation of the air temperature as well [27]. A resolution of the wind speed measurement of 0.01 m s−1 is achievable and a resolution of the wind direction measurement of 0.1° is achievable. The precision and measuring frequency make these anemometers particularly suitable for investigating various turbulence scales. Furthermore, the sensors do not contain any moving parts and are therefore not prone to fouling by airborne dust or sediment.

Fig. 2
figure 2

3D ultrasonic anemometer installed on Queensferry Crossing Bridge [28]

Full size image

In this study, data from five anemometers are used. One of the anemometers is located at top of the central tower (Point A) and the rest are located on both sides of the deck level (Point B and Point C) as shown in Fig. 3. The positions of deck level sensors are all 162.5 m away from the central tower, Point B is on the north of the central tower and Point C is on the south side. On each side of the deck, the west side of Point B will be noted as Point BW and the east side as Point BE (similarly for point C).

Fig. 3
figure 3

Schematic of the bridge and locations of the anemometers used in this study [28]

Full size image

2.2 Data collection and restoration

Raw data available for this research were collected between July 2020 and December 2020. The data were provided to the authors by Transport Scotland in compressed pgdump files. The average size of the compressed raw data for each month was 80 GB. The restoration of the data and its subsequent storage was conducted on the Sonic, the High-Performance Computing Cluster in University College Dublin. PostgreSQL, which is an open-source relational database management system, was used to restore and warehouse the data. The process of schema creation and data restoration was performed in parallel with 20 CPU cores. It took approximately 12 h to complete data restoration for 1 month. Overall, the 6 months of restored data comprise of approximately 13 TB.

2.3 Wind scenarios

The monthly wind data collected by the tower top anemometer (Location A in Fig. 3) are shown on rose plots in Fig. 4. It is apparent from Fig. 4 that the prevailing winds are from the southwest. A selection of specific wind events was chosen (see Table 1), each 1 h in duration, to use as validation studies in this paper. August 2020 and November 2020 were chosen as representative months given that August is the only month where the winds are predominantly from the northeast and November is the month with the highest average wind speeds. The specific times were chosen due to relative stability of both wind direction and wind speed in duration.

Fig. 4
figure 4

Distribution of wind conditions for each month from tower top anemometer data (Location A in Fig. 3)

Full size image
Table 1 Summary of selected dates and wind events for anemometer at Location A
Full size table

A representative sample of data from the top tower anemometers (location A in Fig. 3) can be seen in Fig. 5. Data were recorded at one second intervals and thus the total number of velocity data points in an hour was 3600. Due to the random nature and volatility of natural wind, it was pragmatic instead to consider the average wind speed, which tends to be relatively steadier and more practical for further comparative studies with CFD. First, the 1-min average was determined from the raw data. This was then compared to the 5-min average which is broadly similar and retains most of the characteristics of the 1-min average. Therefore, the 5-min average data are primarily used for comparisons in this study.

Fig. 5
figure 5

Wind velocity data at each record point from selected time on November 15th at a Point A; b Point BW; c Point BE; d Point CW; and e Point CE

Full size image

3 CFD modelling

The details of the 3D CFD model that was developed of the Queensferry Crossing Bridge are presented in this section.

3.1 Geometry

The central tower and central span of the Queensferry Crossing were modelled rather than the whole bridge structure due to limitations of computational resources. The geometry was in full scale with the same dimensions as the real bridge and all secondary structures and cables were included in the model. The bridge tower is 210 m high. Cables on both sides of the central tower hold the deck (see cross section in Fig. 6) which is 788 m in length and 41.6 m in width. Based on the chosen wind events (Table 1), six different model geometries for the six different wind directions were developed by altering the yaw angle of the bridge model. Figure 7a–c shows the geometries for the November wind events (Case A, B and C) where incoming flow is blowing from the right to simulate southwest wind conditions. Similarly, Fig. 7d, e, f shows the geometries for the August wind events (Case D, E and F) where the northeast winds are simulated with flow blowing from the left.

Fig. 6
figure 6

Queensferry Crossing Bridge deck cross section (dimensions in m)

Full size image
Fig. 7
figure 7

Geometries of the Queensferry Crossing for various wind directions: a Case A (200.41 degrees); b Case B (229.89 degrees); c Case C (268.85 degrees); d Case D (75.27 degrees); e Case E (80.55 degrees); and f Case F (85.05 degrees)

Full size image

3.2 Computational domain

Three types of domain boundary region were defined: the inlet face, the outlet face and no-slip walls. According to Liu et al. [29] the boundary in the horizontal flow directions should be a distance of at least five times of structure height (5H) to allow for the establishment of a realistic flow behind the wake region and the buildings included in the computational domain should not exceed the recommended blockage ratio (3%). Due to the slenderness of bridges, the tower height was not the determinant factor to the blockage ratio, so the bridge width was deemed as the primary length in this study. Thus, the inlet face and the outlet face were set at distances of five times and ten times that of the bridge width away from bridge model, respectively, to achieve the blockage ratio less than 1% and to minimize possible errors due to domain size. The size of the computational domain for each geometry was different to optimise the usage of computational resources and the domain size was decided based on the previously mentioned principle. The inlet and outlet faces were 200 m and 400 m away from the bridge, respectively. The front and back faces were both more than 100 m away from the bridge. The height of the domain was three times the bridge height to ensure that the vertical wind profile was able to fully develop. For example, the domains for the November events are shown in Fig. 8. The green face is the inlet face, the red one is the outlet, and the rest of the faces are no-slip walls. The equivalent length and width of the geometry changes after altering the angle of the model; thus, the length and width of the domain changes accordingly. In Fig. 8a, the equivalent length (z direction) and equivalent width (x direction) of the bridge are around 275 m and 730 m respectively, so the length and width of domain are set as 525 m and 1330 m, respectively. After rotating the geometry, as shown in Fig. 8b, c, the equivalent length of geometry increases, and the equivalent width decreases so that the domain size is adapted accordingly.

Fig. 8
figure 8

Domains for November events’ geometries: a Case A (200.41 degrees); b Case B (229.89 degrees); and c Case C (268.85 degrees)

Full size image

3.3 Mesh scheme

The computational mesh scheme was created using OpenFOAM [30] utilities called blockMesh and snappyHexMesh. The BlockMesh utility provides a structured background grid and snappyHexMesh utility works on further defined regional refinement and near-wall region mesh development. Firstly, the background grid composes around 15,000 cubic cells. Each side of a cell is 35 m in length. Further refinement was applied to the nearby region of the bridge geometry. Seven levels of refined cells were stratified, and each level of cell size was half of the previous one. Three buffer layers were assigned between two adjacent levels to smooth the mesh scheme. An overview of the mesh scheme for Case A is illustrated in Fig. 9, which contains a total of 84,705,350 cells. Each level of cells was manually set based on the distance from the bridge model. The cells at Level 7, which are cubes of size 0.273 m, were placed in the region within 1 m of the bridge. Within the 20 m region, Level 6 cells (0.547 m side length) were arranged. Level 5 cells (1.094 m) and Level 4 cells (2.188 m) were 36 m and 56 m away from the bridge, respectively. This cell arrangement resulted in an average y+ value of approximately 150 being achieved at the wall of the bridge geometry, which is within the optimal range of y+ values (30 to 300) when using the wall function. According to Han et al. [31], the mesh scheme with minimum mesh size of 1 m and an average y+ value of 300 was sufficient for a large full-scale terrain wind fields study. Therefore, the finer mesh scheme was deemed to satisfy the simulation requirement. A mesh sensitivity study was conducted, the results of which will be discussed in Sect. 4.1.

Fig. 9
figure 9

Slice view of mesh: a general view of the mesh grid; and b enlarged view of the mesh grid

Full size image

3.4 WRF model

To determine realistic inlet conditions for the CFD model, the numerical atmospheric system was simulated using a Weather Research and Forecasting (WRF) model, which is a non-hydrostatic, compressible model. WRF models make real-time numerical weather predictions, can model weather events and atmospheric-process studies, perform regional climate simulation, atmosphere–ocean coupling and idealized atmosphere studies [32]. The extent of the horizontal domain of the WRF model is illustrated in Fig. 10. Horizontal grid points are spaced 12 km apart, with 100 grid points in both the east–west and north–south directions. There are 60 vertical levels, with the highest around 20 km above sea level. The actual location of the Queensferry Crossing Bridge and its nearest weather station (Edinburgh airport) along with the closest grid point from which WRF output was recorded are illustrated in Fig. 10b. The time resolution of the WRF model output was in a second. Data from Edinburgh Airport (EGPH) weather forecast station (the actual location is highlighted as a red star in Fig. 10b) were used as reference to validate and adjust the WRF model by comparing simulated results at the closest grid point to the EGPH (red circle in Fig. 10b) with official records (available open source from the UK Automated Surface Observing System [33]). Then, wind velocity data from selected dates at the closest grid point to the bridge location (black star in Fig. 10b) were collected and processed as initial velocity condition for CFD simulations (see Fig. 11). Given that the bridge location is in an open sea area, the spatial distance effects between the locations of the model grid point and the bridge can be negligible. The turbulence kinetic energy (k) and specific turbulence dissipation rate (ω) were calculated according to velocities from sea level to corresponding height of the simulation.

Fig. 10
figure 10

Simulation domain for the WRF model: a overall WRF domain; and b closest grid points from which WRF output was recorded

Full size image
Fig. 11
figure 11

Incoming wind profile for Case A obtained from WRF model

Full size image

3.5 CFD governing equations

The main governing equations for the fluid flow in the CFD analysis are the mass conservation and momentum equation, which together constitute the Navier–Stokes equations. However, solving the time-dependent Navier–Stokes equations for fully turbulent flows is only feasible for very simple cases due to limitations in computing power. Meanwhile, for the complicated cases of real-world structures, it would be computationally expensive to accurately simulate and describe every eddy in a flow. Therefore, time-average transport equations, called Reynolds-averaged Navier–Stokes (RANS) equations, are widely applied. In this study, the simulative flow is set to be incompressible, isothermal, Newtonian and statistically steady. Gravity effects are neglected. Equation 1 is the mass conservation equation and Eq. 2 is the momentum equation:

$$\nabla \cdot U = 0 ,$$
(1)
$$\nabla \cdot \left( {UU} \right) = \nabla \cdot \left[ {v_{{{\text{eff}}}} \left( {\nabla U + \left( {\nabla U} \right)^{{\text{T}}} } \right)} \right] - \nabla p,$$
(2)

where U is the time-averaged velocity, p is the time-averaged kinematic pressure, and \(v_{{{\text{eff}}}}\) is the effective kinematic viscosity, which is the summation of laminar kinematic viscosity (1.47 × 10−5 m2/s) and the turbulent kinematic viscosity (calculated with k − ω SST turbulence model).

OpenFOAM [30] uses the finite volume method to discretise terms in the RANS equations and solve these partial differential equations. The equations are integrated over the control volumes (cells) of the mesh. The volume integral of the divergence are converted to surface integrals by using Gauss’s theorem. Then, by summing the value from every face of a cell, the surface integrals of a vector field can be obtained.

3.6 Boundary conditions

On the inlet patch, parameter values for U, k and ω were mapped from the WRF model onto at least 4424 points based on the size of the inlet patch. The distances between each point are vertically 2 m and horizontally 35 m to precisely prescribe inlet boundary condition as the WRF simulation results. The turbulence intensity can be obtained with the WRF model results using following equation [34]:

$$I = u^{\prime}/U ,$$
(3)

where μ′ is the root-mean-square of the velocity fluctuations. The inlet–outlet boundary condition is applied at the outlet patch to prevent reverse flow. The remainder of the wall regions was set to be no-slip walls and adaptive wall functions were applied for k and ω configurations. The process strategies can automatically switch between high and low Reynolds number based on the y+ value. A summary of boundary conditions is shown in Table 2.

Table 2 Summary of boundary conditions of CFD simulations
Full size table

3.7 Numerical configuration

Simulations in this study employed the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm [35] to solve the pressure–velocity coupling. Gradient and Laplacian terms were discretized using Gaussian integration with linear interpolation. Convection and advection terms were discretized using a second-order accurate linear-upwind scheme. The momentum and pressure equations were solved in an iterative manner. First, an incipient velocity field can be obtained as the momentum equation is initially solved. Then, the pressure equation can be solved so that the pressure field is modified accordingly. After that, the velocity field was calculated and updated again based on the latest pressure field. K and ω were also calculated to update the turbulent viscosity. Thereafter, fluxes at the cell faces and the boundary conditions were updated accordingly, and the above process was repeated until convergence was obtained. Drastic oscillation of the residuals can happen, which can lead to corruption of simulations. Therefore, relaxation factors were set as 0.4, 0.5 and 0.7 for P, k/ω and U, respectively, to balance the stability and process speed. The CFD model was considered converged at 1000–1500 iterations based on the residual value remaining constant around 10−4. The simulations were allowed to run for up to 2000 iterations to confirm that the residuals remained constant beyond 1500 iterations.

4 Validation of the CFD models

4.1 Mesh sensitivity study

The quality of a mesh is one of the primary factors that can severely affect the accuracy of CFD simulation results. Theoretically, a perfect mesh scheme could be achieved with unlimited computational resources. However, due to the inevitable limit of computational power, it is meaningful to investigate an efficient level of mesh refinement to obtain relatively accurate results with practical mesh schemes. Several existing full-scale numerical works on terrain studies [31, 36] relating to wind field data had simulated velocity profiles with average differences around 30% compared to field measurement, and the maximum average differences could reach 45%. Therefore, in this study, good accuracy was deemed as (1) the average difference on each point is no more than 20% in each case; and (2) the maximum difference is less than 35%.

4.1.1 Effect of changes in refinement level

First, three levels of refinement are compared for Case A, Case B and Case C. As shown in Fig. 9, the highest level of cells near the bridge model was Level 7 (0.237 m side length). In this study, 6 levels and 8 levels of refinements are structured in a similar way as that described in Sect. 3.3. For the refinement consisting of 6 levels, the highest level of cells near the bridge model was Level 6 (0.547 m side length). For the 8 levels of refinement, the near bridge region was surrounded by Level 8 cells (0.118 m side length). All the other configurations such as domain size, boundary conditions, etc., remain same for each simulation.

In order to create an equivalent environment with field monitoring, every numerical timestep is deemed as a second. The average value of every sixty timesteps is regarded conceptually equal to 1 min average field data. Thus, averaged values of every five averaged simulation results were used to compare with five minutes averaged field data.

As illustrated in Fig. 12, results of wind velocity for Case A on three different refinement levels are plotted and verified with processed field data (referred to 5-min average value in Fig. 5). Results from the different refinement levels only show a slight difference (as shaded bands in Fig. 12). Generally, the difference among CFD results from different refinement levels is less than 10% for Case A. The simulation results show a similar trend with field data as well, especially for the deck level results (Fig. 12b–e). In Fig. 12c, e, all three groups of CFD results have an average relative difference of less than 15%. The maximum average difference compared to field data was 12.3% which was from the Level 6 refinement results (the average velocity difference was 0.66 m/s). As for the west side records (Fig. 12b, d), the maximum average difference was from the Level 6 scheme as well, which is 23.6% (1.89 m/s). The Level 8 mesh scheme provides the closest wind velocity at each time point to the field data, which suggests that finer refinement tends to obtain lower differences when compared to the field data on deck level. In Fig. 12a, it is interesting to note that the Level 8 scheme is overestimating the wind velocity on the top of the tower. The average difference was 11.6% (1.64 m/s) for the Level 8 scheme, which was higher than the other two schemes. The Level 6 scheme obtains the closest results when compared to field data as the average difference was 10.3% (1.45 m/s). These simulation results showed agreement with conclusions suggested by Van Hooff et al. [37] that RANS turbulence models tend to underpredict the turbulent kinetic energy in areas above and below the incoming flows. As for the k − ω SST turbulence model, it tends to provide lower levels of turbulent kinetic energy in the incoming stagnation region and larger velocity gradients in the shear layer of flows compare to the experiments. Slice views of wind velocity contour on near bridge region for Case A in level 7 mesh scheme are illustrated in Fig. 13.

Fig. 12
figure 12

Effects of changes in mesh refinement level for Case A: a velocities at Point A; b velocities at Point BW; c velocities at Point BE; d velocities at Point CW; and e velocities at Point CE

Full size image
Fig. 13
figure 13

Wind velocity contour on near bridge region for Case A in level 7 mesh scheme: a slice view at Point A; b slice view at Point B and c slice view at Point C

Full size image

Table 3 illustrates the results comparison for Case A, Case B and Case C, respectively. CFD results were able to show good agreement with field data. For Case B, results show similar features to the results from Case A. Level 8 scheme has the closest estimation of the wind velocity for the tower top level, but velocity on the deck level is overestimated especially for the west-side results. One of the potential reasons could be due to the traffic influence on the deck level since the time of Case B occurred during a weekday at 6.00–7.00 h rather than in Case A, which occurred at the weekend from 20.00 to 21.00 h. As per the statistics report from Transport Scotland [38], the average traffic flow on M90 in a weekday was 21.7% greater than at the weekends. The average traffic flow on M90 in early morning was around 35% higher than at night-time in weekends. Therefore, the traffic effect which was neglected in the simulations could be a factor on the prediction accuracy of the deck level. For Case C, results from the Level 7 scheme and the Level 8 scheme can provide closer prediction than the Level 6 scheme. Overall, CFD results demonstrate that numerical simulations can perform sufficient estimation of the wind velocity, and at some locations, even similar fluctuations can be simulated when compared with field data. As expected, higher levels of mesh refinement can provide higher resolution outcomes; however, higher resolution may trigger overestimate issues which can lead to larger difference between CFD results and field data.

Table 3 Summary of the mesh refinements study: results difference between field data and CFD results for Case A, B and C
Full size table

Table 4 summarizes the comparison of field data and CFD results for Case D, E and F using only Level 7 mesh refinement. The maximum difference, minimum difference and average difference of each monitoring point are shown. Case F performs best as the average difference of every monitoring point is all below 10%. Performances from the other two cases are good as well. It is worth noting that the average difference is no more than 20% in each case, meanwhile, the maximum differences is less than 35% suggests that the fluctuations of wind velocity are reasonably well captured by CFD.

Table 4 Summary of the comparison between field data and CFD results for Case D, E and F using Level 7 mesh refinement
Full size table

Table 5 summarizes the total number of cells in each of the mesh schemes. The Level 6 mesh schemes contain around 55 million cells, the Level 7 mesh schemes contain around 85 million cells and the Level 8 schemes around 120 million cells. The average subsequent computational time for all the Case A, B and C at each mesh Level is presented in Table 6. It can be seen that the clock time of simulations with the Level 8 mesh is more than 1.5 times of that using the Level 7 mesh. The Level 8 mesh schemes are able to obtain 1.9% closer prediction to the field data on average, however, the possibility of overestimation and the significant cost of computational power makes the Level 8 scheme less efficient than the Level 7 schemes. Therefore, the Level 7 mesh scheme was chosen for the remainder of this study as it has the ability to perform relatively accurate prediction results and it is relatively more efficient.

Table 5 Overall cell number of the three mesh schemes
Full size table
Table 6 Average parallel simulation time of simulations for Case A, B and C
Full size table

4.1.2 Effect of changes in refinement distance

Using the Level 7 mesh scheme, a refinement distance study was conducted to investigate the trade-off between savings in computational power (fewer high level cell regions) and accuracy of results. In this investigation on refinement distance, all analyses have 7 levels, but within those levels refinements in the mesh were made in order to reduce computational time, i.e., reduce the time from 21.7 h in Table 6. This study was carried out on Case A and all the numerical configurations are the same except for the refinement distance. The original mesh scheme distance arrangement is illustrated in Fig. 8b, another two different refinement distance mesh schemes were created, and details of refinement distance can be found in Table 7. For Mesh No. 1, the distance of each level cells (except level 7 cells) to the geometry surface is half of the original scheme, and Mesh No. 2 further reduces the distance to the bridge surface of each level. In general, narrowing the distance of high-level refinements can significantly save computational resources, but as expected, at the expense of accuracy of simulation results to a certain extent. There appears to be a minimal effect at the top of the tower to changes in refinement distance. However, when it comes to the deck level, the capability of the lower level mesh at describing the change of flow is not sufficient when compared with the field data. The maximum difference is greater than 35% for Mesh No. 1. Mesh No. 2 shows much higher average results greater than the target criteria (20%). It is worth mentioning that both Mesh No. 1 and No. 2 had average y+ values of bridge geometry larger than 200, which is higher than the original mesh (150). The minimum y+ value of the bridge in the original mesh was 1.1, which was lower than Mesh No. 1 and No. 2. These values suggest that the original mesh scheme is sufficient and should be applied for the remainder of the study. Moreover, based on the mesh density results, a minimum y+ value around 1 and an average y+ value around 150 can be recommended for full-scale bridge aerodynamic CFD studies.

Table 7 Summary of refinement distance study for Case A
Full size table

4.2 Domain study

It is typical for CFD domains to be created in the shape of rectangles. One of the reasons may be that the results of the majority of CFD studies on bridges are compared with results from wind tunnel studies which tend to be rectangular in shape. Another reason may be that, in rectangle domains, it is easier to control the size of the background grid and side length of cells. However, rectangular domain shapes are inconvenient for some yaw angles studies. One of the most typical ways to achieve different yaw angles is by rotating geometries [39, 40], which is a similar method to angle of attack studies. But changes in yaw angles of a model are likely to face issues such as unsuitable domain size/refinement region size and changes of desired data point locations. For these reasons it is common to have to recreate a domain, a mesh scheme and adjusted other parameter settings accordingly. Another common approach is to apply wind flows from two directions. Different yaw angle wind can be achieved by adjusting the relative velocity of the two flows [41, 42]. Although the flexibility and logic of this method are widely accepted [41], the potential turbulence motion caused by two flows should not be ignored. Therefore, in this study, a cylinder-like domain is designed and investigated. The results are then compared with more common rectangle domain results under similar configurations.

As can be seen in Fig. 14, instead of rotating the bridge geometry to create different yaw angle cases, in a cylindrical domain, different yaw angles can be easily achieved by changing positions of the inlet face (green face). The radius of the domain is 500 m to ensure that there is enough space between the geometry and outlet in the majority of angles for the development of turbulence. The height of the domain is the same as the typical rectangle domain, which is 630 m. The inlet face is a plane rather than an arc face to avoid interlacement of the incoming flow. The size of the inlet face was based on the equivalent length of the bridge under a certain yaw angle. Moreover, mesh schemes for the cylindrical domain case were mostly the same as the rectangular domain cases. Normally, o-grid blocking is a common mesh method for cylindrical domains when fluid direction normal to the radial plane [43]. However, due to the long shape of the bridge, non-uniform and non-orthogonal structured grids are inevitable in o-grid mesh generation. According to Volk et al. [44], o-grid refinement on the cylindrical domain produced a rapid increasing in the number of grids within the desired range, but the results over the grids error rate showed large oscillation. A cubic or hexahedron mesh scheme tended to result in a significant decrease in error rate compared to a non-uniformed o-grid mesh scheme. Thus, the cubic-based mesh scheme was utilised in this study and most of cells are intended to be set with 35 m side lengths (Fig. 15a). The refinement distance was enlarged based on size of the domain. Level 7 cells were kept within the 1 m region near the bridge model, then, Level 6 cells were arranged until 30 m away from the bridge. Level 5 and Level 4 refinements were 60 m and 90 m away from the bridge, respectively (Fig. 15b).

Fig. 14
figure 14

Cylindrical domains for November events: a Case A (200.41 degrees); b Case B (229.89 degrees); and c Case C (268.85 degrees)

Full size image
Fig. 15
figure 15

Cylindrical domain: a overview of background grid; and b slice cut view of the mesh scheme

Full size image

The numerical configurations were the same as the previous simulations in Sect. 3. Simulation results of Case A, B and C are compared with the rectangular domain results. In Case A (Fig. 16), results determined by the cylindrical domain are closer to field data than the rectangular domain results. Especially on the tower top level, the cylindrical domain simulation produces on average a 7.7% better prediction than the rectangle domain study, which only had 3.4% difference on average (0.48 m/s). At deck level, results from both domains are similar. Within each simulation, there are times when the cylindrical domain gives closer results to the field data, and other times when the rectangular domain does. However, in Case B and Case C, the results of the cylindrical domain model are not as close to the filed data compared with those from Case A. As illustrated in Table 8, most wind velocity determined by the cylindrical domain simulations have moderately larger differences than the rectangular domain. The maximum difference is as high as 49.1% which happens at Point CE in Case C. Besides, it is also worth noting that at some recorded points, the cylindrical domain simulation is able to obtain results at the same level of accuracy as the rectangular domain simulations. For instance, in Case B, results from the cylindrical domain are closer to the field data. The average difference to the field data at deck level is all lower than 10%. However, in Case C, results from the cylindrical domain model suggests that the performance of the model can be affected under a certain extent changes of the yaw angle.

Fig. 16
figure 16

Domain study: comparison of wind velocity at each record point for Case A

Full size image
Table 8 Summary of the domain study: results difference between field data and CFD results for Case A, B and C using Level 7 mesh refinement
Full size table

In general, this study suggests that the cylindrical domain could potentially be a good substitution for rectangular domains in studies involving various yaw wind angles. The cylindrical domain has the ability to provide reasonable numerical estimation of wind velocities. Nevertheless, the robustness of the cylindrical domain model still needs to be further investigated. Due to the curved shape of a cylinder, the mesh grid may have some irregular shapes which can have a negative influence on numerical calculations, or it could affect the mesh quality potentially causing floating point problems. Therefore, further study and improvements are necessary before widespread implementation of cylindrical domains.

5 Conclusions

In this study, for the first time, full-scale bridge aerodynamic CFD simulations were developed with the k − ω SST turbulence model and compared to field data. A 3D full-scale model of the central tower section of the Queensferry Crossing Bridge, United Kingdom, was created and used in this study. To replicate the real-world wind conditions on the site of the Queensferry Crossing, a full-scale WRF model provides the initial wind velocity conditions of the selected time and dates for the CFD simulations to build atmospheric boundary layer inflows. Then, simulation results were validated with field monitoring data namely accelerometers mounted on the tower and deck. The numerical predictions were able to describe fluctuating wind profiles with relative differences of less than 20% to the field data in 1-h durations, which indicated that the model accurately captures 3D effects and the full scale WRF simulation provides good estimates for input wind velocity conditions. Meanwhile, a mesh density study was conducted to investigate simulation performance under three different refinement levels, which brought smallest cell size to 0.547 m, 0.237 m and 0.118 m, respectively. The effect of three different refinement areas was then investigated and a mesh scheme with appropriate refinement distance was used for further studies as it showed the optimum balance of efficiency and accuracy. Based on this mesh density analysis, a minimum y+ value around 1 and an average y+ value around 150 can be recommended for future full-scale bridge aerodynamic CFD studies.

Lastly, a cylindrical-like domain was designed and applied for this first time to a full-scale bridge model. The cylindrical-like domain aimed to provide a convenient and reliable way to investigate bridge aerodynamics regarding changes in yaw angle. The results from the cylindrical domain analyses were compared with traditional rectangular domain analyses results and showed good reliability.