Analyzing Wind Effects on Long-Span Bridges: A Viable Numerical Modelling Methodology Using OpenFOAM for Industrial Applications

Abstract

Aerodynamic performance is of critical importance to the design of long-span bridges. Computational fluid dynamics (CFD) modelling offers bridge designers an opportunity to investigate aerodynamic performance for long-span bridges during the design phase as well as during operation of the bridge. It offers distinct advantages when compared with the current standard practice of wind tunnel testing, which can have several limitations. The proposed revisions to the Eurocodes offer CFD as a methodology for wind analysis of bridges. Practicing engineers have long sought a computationally affordable, viable, and robust framework for industrial applications of using CFD to examine wind effects on long-span bridges. To address this gap in the literature and guidance, this paper explicitly presents a framework and demonstrates a workflow of analyzing wind effects on long-span bridges using open-source software, namely FreeCAD, OpenFOAM, and ParaView. Example cases are presented, and detailed configurations and general guidance are discussed during each step. A summary is provided of the validation of this methodology with field data collected from the structural health monitoring (SHM) systems of two long-span bridges.

1. Introduction

Assessment of the wind effects on long-span bridges is a critical aspect of their safe design. When wind blows on a bridge, it exerts a distributed pressure force on the structure, which is related to the velocity of the wind. Forces are resisted by bridges and transferred to the ground through mechanisms of bending, torsion, and shear. In exposed sites, wind loads can sometimes have a greater impact than any other bridge load, such as thermal loading and traffic loading [1]. An essential role of any bridge engineer is to estimate the effects of prevailing wind conditions on a bridge and suitably design the various bridge elements for this loading, yet the forces that bridges experience in high winds are poorly understood. More specifically, when designing long-span bridges, engineers can face aerodynamic and aeroelastic challenges due to wind effects, such as buffeting [2], flutter [3], vortex-induced vibrations [4], galloping [5], and cable vibrations [6]. To predict and estimate the influences of these problems on bridge design, several key parameters are of interest, namely the aerodynamic coefficients, vortex shedding frequencies, and flutter derivatives. The aerodynamic coefficients are typically the most frequently used. They are dimensionless coefficients that are used in calculations of aerodynamic forces and flutter derivatives.

The accepted state of the art is that bridge design is based on wind tunnel testing [7,8,9,10]. Although wind tunnel testing for design offers data, the data are limited due to assumptions of stationarity and scalability, particularly of Reynolds and Strouhal numbers. Furthermore, locations of sensors must be decided a priori and can only provide information at discrete measurement points. Flow visualization of the entire domain (e.g., PIV) is typically prohibitively complicated and expensive. Scheduling time at a wind tunnel facility can prove expensive both in terms of time and cost, particularly when repeat tests might need to be conducted. Aerodynamic admittance, based on the seminal work of both Davenport and Scanlan, is assumed to act strip wise (1D) for bridges to create the buffeting forces, and only recently have new analytical forms been proposed [11] to reflect the growing length and complexity of bridges. This increasingly means that 1D and zero-yaw angle winds are no longer reasonable assumptions, and together with the previous limitations, this makes wind tunnel testing hugely limiting for this application.

The opportunities afforded by the advent of greater computational capacity offers the use of computational fluid dynamics (CFD) as an alternative method [12,13,14,15,16,17]. Even in the early stages of development, when computational power was significantly limited, some pioneers [18,19,20,21,22,23] demonstrated the capacity of CFD simulations in this area. CFD has been used to assess the non-uniform wind characteristics in mountainous areas and its effects on the flutter performance of a long-span suspension bridge [24]. In July 2008, a benchmark study on the aerodynamics of rectangular cylinders was initiated with the aim of analyzing the turbulence and separated flow around an elongated rectangular cylinder, analogous to a bridge [25,26,27,28,29,30]. Analyses have used both CFD and WTT, and results have yielded valuable insights. Computational methods have been used to investigate the edge degree-of-sharpness and free-stream turbulence scale effects on the aerodynamics of a bridge deck [31] to investigate the impact of gap distance between twin-box decks [32], as well as to optimize the shape of bridge decks considering aeroelastic and structural constraints [33,34]. More recently, authors have begun to investigate the effects of bridge deck movement as part of these computational analyses [35,36,37,38].

Currently, the Eurocode on wind actions [39] provides guidance on calculating wind effects under specific conditions with sets of equations and recommends wind tunnel testing in all other conditions. CFD is not yet specified in the Eurocode as a means of estimating wind effects on long-span bridges. Interestingly, the Eurocodes are currently undergoing revision, and the latest draft includes an annex permitting the use of CFD in certain circumstances [40]. However, there is a lack of prescriptive guidance to accompany these design codes for bridge design engineers. There is demand from bridge designers, bridge owners, and bridge operators for a CFD methodology that they can use in assessing their existing bridge stock for adverse impacts due to prevailing wind conditions, as well as assessing newly planned bridges [41]. According to Zhang et al. [42], many publications that present analysis of wind effects on bridges using CFD have provided descriptions of their methodology in relatively little detail and rather choose to focus on analyses of the results [43,44,45]. Most currently available online guidance does not consider the unique application of wind effects on long-span bridges. The National Standards Authority of Ireland is currently leading the Irish participation in revisions to the Eurocode on Wind, and therefore this paper directly supports that work. It explicitly articulates an approach to CFD model development that is sufficiently generic as to accommodate many different forms of long-span bridges and that is viable for industrial applications, with enough detail for bridge designers to use in studying wind effects on a sample long-span bridge using OpenFOAM. Crucially, this methodology has been tested and validated with structural health monitoring (SHM) data collected from two long-span bridges in the author’s previous work on the Queensferry Crossing in Edinburgh [22,23] and the Rose Fitzgerald Kennedy Bridge in Ireland [21,46]—further details provided in Section 11. While more sophisticated and more computationally demanding modelling (such as a large eddy simulation incorporating fluid structure interaction) might yield more accurate results, this methodology has been demonstrated to be sufficiently accurate when compared to field data, with a level of computational capacity that makes it viable in industrial applications.

The workflow presented in this paper starts with defining some modelling assumptions in Section 2 followed by developing the geometry, as described in Section 3. After sketching and exporting the 3D geometries to STL files, the configuration of an OpenFOAM case is initiated by creating the computational domain and the mesh (described in Section 4). Next, in Section 5, boundary conditions are specified for different patches. The paper then presents the configurations for multiple dictionaries related to discretization (Section 6), the algorithm solvers (Section 7), and parallel decomposition (Section 8). The process of conducting a mesh sensitivity study is then presented in Section 9. Finally, Section 10 describes the post-processing steps where both quantitative and visualized results are considered. Figure 1 provides an overview of the workflow proposed in this paper.

2. Modelling Assumptions

In an analysis of a typical bridge, the wind velocity at the site will normally have a Mach number that is much smaller than 0.3, and so an assumption of incompressible flow is valid for CFD analysis [47]. Furthermore, the temperature and gravity effects in bridge aerodynamic studies are often neglected, as their contribution to the forces are negligible. During bridge design, aerodynamic coefficients are often used to estimate the aerodynamic performance of a bridge [48,49,50]. These parameters can be calculated from simulations using the Reynolds-averaged Navier–Stokes (RANS) formulations (for example, as implemented in simpleFoam):

∇ · U = 0,

(1)

∂U/∂t + ∇ ⋅ (UU) = ∇ ⋅ [ν_eff (∇U + (∇U)^T)] − ∇p

(2)

where U is the Reynolds-averaged velocity, p is the Reynolds-averaged kinematic pressure (pressure divided by density), and ν_eff is the effective kinematic viscosity, which is the summation of laminar kinematic viscosity (ν) and the turbulent kinematic viscosity (ν_t). The value of the laminar kinematic viscosity, ν, can be configured in the transportProperties file in the constant/directory, which is related to the temperature of the flow. For instance, ν has a value of 1.47 × 10⁻⁵ m²/s at a temperature of 15 °C. In Equation (2), the ∂U/∂t term accounts for transience in the system (URANS) but is omitted for steady-state RANS simulations.

The unknown ν_t in RANS and URANS simulations can be calculated using various turbulence models. The use of the k-ω SST turbulence model [51] has become popular in the study of bridge aerodynamics [52,53,54] due to its ability to model turbulence behavior in the boundary layer [55]. It should be noted that there are more complex models that might deliver more accurate results, such as the four-equation $\overline{v^2}$-f model [56], the four-equation γ-Re_θ model [57], and the seven-equation SSG/LRR-ω full Reynolds stress model [58]. Since these models contain more transport equations than the two-equation k-ω SST model, they often require significantly more computational power and are more challenging to converge. Furthermore, the implementation of some of these models, such as the γ-Re_θ, also needs empirical correlations to configure inlet parameters, which requires greater effort in calibrating the model to fit specific engineering problems. Given that the purpose of the workflow in this paper is for the engineering design of long-span bridges, computational affordability is of significant importance. In this application, the k-ω SST turbulence model is both affordable and reliable and hence has been proposed for the application of modelling wind effects on bridges in RANS and URANS simulations. To solve the aforementioned governing equations with the k-ω SST turbulence model in OpenFOAM, the simpleFoam solver should be used for RANS simulations, and the pimpleFoam or pisoFoam solver should be used for URANS simulations. In the example case of URANS simulations, the pimpleFoam solver is used because it allows for large time steps, which potentially speeds up the transient simulations. It should be noted that if accurate predictions of turbulence structures are of interest, large-eddy simulations (LESs) [59] and detached-eddy simulations (DESs) [60] are recommended. These simulations can also be performed using the pimpleFoam or pisoFoam solver. However, compared to URANS, LESs and DESs require significantly more computational power. Figure 2 visually depicts key elements of the framework, including boundary conditions, the turbulence model, and discretization methods.

3. Geometry Development

A geometry file describing the bridge structure is a required input of OpenFOAM. There are several software options that can generate 3D bridge geometries. One option is to generate 3D geometries based on 2D sketches using AutoCAD version 2021, commonly used software in civil engineering designs. Based on the authors’ experience, 3D solids sketched in AutoCAD often contain redundant lines and points and sometimes unclosed curves, which affect the quality of the geometry. The graphic user interface (GUI) of AutoCAD was designed to work in 2D and is not user-friendly when sketching in 3D. The authors recommend using software that adopts a parametric design approach, such as SOLIDSWORKS [61] and FreeCAD [62]. While SOLIDSWORKS 2021 version 29 offers the most mature and reliable solutions for industrial design, it is a commercial software that requires a license fee. In contrast, the Python-based software, FreeCAD, is an open-sourced alternative that offers the parametric design approach and a large suite of highly customizable modules. Figure 3 shows a 3D geometry of an example case generated by FreeCAD (version 0.19). The model contains five components, namely vertical handrail bars, horizontal handrail bars, fascia beams, central barriers, and the bridge deck. These components were sketched separately and then assembled. The bridge deck geometry was then configured with different pitching angles (from −10 degrees to +10 degrees) to prepare for simulations of the bridge decks under different angles of attack (Figure 4).

After creating the geometry, users need to export the file into a format readable by OpenFOAM. While OpenFOAM can read geometry files of various formats, the authors recommend exporting geometries as ASCII STL files, where keywords within the files can be edited by text editors.

4. Creating the Mesh

There are various approaches to creating the internal mesh for computing in OpenFOAM. For example, structured meshes can be generated using the blockMesh utility. Unstructured meshes can be generated by the snappyHexMesh utility or the cfMesh library. Alternatively, users can import meshes from third-party mesh generators, such as ICEM [63] and Gmsh [64]. However, 3D geometries of the bridge might contain superfluous details, and therefore an unstructured mesh can be more practical. For this application, the bridge model is contained within a hexahedron-dominated unstructured mesh generated from the snappyHexMesh utility of OpenFOAM. Firstly, the blockMesh utility is used to define the computational domain and generate the structured background mesh. Figure 5 shows the dimensions and patch names of the domain for an example case. The domain dimensions presented here are typical of the size of wind tunnels used for bridge assessments. The authors propose that users should select a domain size that is relative for their application so as to maintain the blockage ratio below 5% [65,66].

Upon creation of the domain and background mesh, snappyHexMesh is used to perform cell refinement. The development of a mesh is not straightforward as there are various quality requirements. Firstly, the mesh should be fine enough at the bridge surface so that details of the bridge geometries are resolved and accounted for in the simulation [67]. Secondly, the mesh boundary layer at the bridge surface should be fine enough to provide a mean y⁺ value below 1 so that the use of k-ω SST turbulence is valid [51]. This is especially important for simulations with geometries containing small details. In these cases, cells within the near-wall region are refined to fully resolve the geometric details and so often correspond to a small y⁺ value of less than 30, i.e., low-Reynolds-Number (Low-Re) wall condition [68]. While OpenFOAM provides wall functions that can be adopted in Low-Re wall condition [69], wall functions do not perform well in the buffer zone where y⁺ ranges from 5 to 30 [70]. Therefore, the most commonly adopted approach is to refine the boundary layer so that the buffer zone is avoided and the viscous sublayer is resolved by the mesh, i.e., the y⁺ value is less than 1. Thirdly, the maximum non-orthogonality and the maximum value of skewness of internal and boundary cells within the mesh should be below 65, 4, and 20 degrees, respectively [67]. Finally, simulations with the selected mesh scheme should also be affordable. More specifically, the mesh is expected to have a cell count of less than 64 million, estimated by multiplying the maximum number of CPU cores (128) used for the study with the maximum number of cells proposed to be distributed on each CPU core (500 thousand). It would be expected that RANS simulations of this type of study with a limit on the cell count would be expected to be completed within 20 h for approximately 8000 SIMPLE iterations. However, transient simulations require significantly more computational efforts, and so URANS simulations of 1 s in this study with the same mesh as in RANS simulations are expected to be completed within a week. Due to these requirements, the development of a suitable mesh is an iterative process. Here, Figure 6 shows a proposed final mesh developed for this example case. A mesh sensitivity study for this example is presented later in Section 9.

As shown in Figure 6, cells within the background mesh have a thickness of 0.05 m. In this case, the y+ values at patches of the domain walls are approximately 250 to 300, which is valid for the use of wall functions [67]. The y⁺ value at the bridge surface should be lower than 1. In this example, a refinement of six levels is applied using snappyHexMesh, where five buffer layers were added between adjacent levels to provide a smooth transition between cells of different sizes. In addition, to better capture the flow characteristics in the region that is close to the geometry, a refinement box was used where cells were firstly refined to level 4 and then refined to level 6 at the surface of the geometry. This led to a thickness of approximately 7.8125 × 10⁻⁴ m for the smallest internal cells, so the smallest geometrical detail with a thickness of 1.0 × 10⁻³ m would be captured. Afterwards, the addLayerControl in snappyHexMeshDict was enabled to insert eight layers of boundary cells. All of the above schemes would lead to a mesh with approximately 30 million cells, with most of them hexahedrons. In this case, the y⁺ value at the bridge surface was determined to be 0.45. The aspect ratio was 3.47.

5. Boundary Conditions

Boundary conditions are specified by configuring the dictionary file of each parameter in the 0/directory.

Table 1. Summary of boundary conditions applied in the example case.

Patch Name	Parameter	Type of Condition in OpenFOAM	Value	Unit
Inlet	U	fixedValue	P1	m/s
	p	zeroGradient	-	m2/s
	k	fixedValue	P2	m2/s
	ω	fixedValue	P3	s−1
	νt	calculated	-	m2/s
Outlet	U	inletOutlet	-	m/s
	p	fixedValue	0	m2/s2
	k	inletOutlet	-	m2/s2
	ω	inletOutlet	-	s−1
	νt	calculated	-	m2/s
Side walls	U	noSlip	0	m/s
	p	zeroGradient	-	m2/s2
	k	kqRWallFunction	-	m2/s2
	ω	omegaWallFunction	-	s−1
	νt	nutkWallFunction	-	m2/s
Bridge	U	noSlip	0	m/s
	p	zeroGradient	-	m2/s2
	k	fixedValue	1 × 10−	m2/s2
	ω	omegaWallFunction	-	s−1
	νt	nutLowReWallFunction	-	m2/s

presents a summary of boundary conditions for each parameter at different patches for this example case of modelling wind effects on bridges.

At the inlet patch, fixed values of P₁, P₂, and P₃ are given to U, k, and ω, where P₂ and P₃ are related to P₁, turbulence intensity (I), and turbulence length scale (l) and can be calculated using Equations (3) and (4):

k = 3/2 (UI)²

(3)

ω = 0.09 − 1/4√k/l

(4)

In engineering practices, these values can be determined from multiple sources. One option is that there might be an anemometer installed near the bridge site, from which P₁ can be determined for the period of interest, and I and l can subsequently be determined. Another possibility is that there may be a nearby weather station from which data can be obtained. A third option is to perform complex meteorological simulations based on weather data, the results of which can also provide P₁, I, and l. In some cases, for bridge construction projects, prior wind tunnel tests are sometimes performed, from which these values could be taken.

It should be noted that this workflow omits consideration of the atmospheric boundary layer (ABL). Should the ABL be of concern, say in a full-bridge simulation that includes terrain effects, U, k and ω can be configured with conditions of atmBoundaryLayerInletVelocity, atmBoundaryLayerInletK, and atmBoundaryLayerInletOmega, respectively, at the inlet patch. In this workflow, at the outlet patch, the inletOutlet condition is configured for U, k, and ω to prevent the reverse flow, which could be replaced by the freestream condition if the ABL is to be modelled. At wall patches of the computational domain, corresponding y⁺ values are larger than 30, and hence kqRWallFunction, omegaWallFunction, and the nutkWallFunction are used for k, omega, and nut, respectively. At the patch of the bridge deck, the y⁺ value is below 1, so the viscous sub-layer is resolved within the boundary cells. Theoretically, k can be configured with a value of 0 m²/s². But it is configured with a small value (1 × 10⁻¹⁰ m²/s²) to avoid float point issues related to division by a zero. It is noted that omegaWallFunction in OpenFOAM is a linear wall condition that fits both high-Reynolds-number and low-Reynolds-number conditions, and so it is configured for ω at the bridge patch. In contrast, the nutLowReWallFunction is configured for ν_t at the bridge patch, which, according to Liu [69], is designed for ν_t at the low-Reynolds-number condition.

6. Discretization of the Equations—fvSchemes

All terms in the governing equations are discretized using the nominally second-order cell-centered finite volume method [71] except for terms in turbulence models. Figure 7 shows the content of the fvSchemes file in the system directory of each RANS simulation case, where detailed numerical schemes of each term of the governing equations are configured. Gradient terms are discretized using Gaussian integration with linear interpolation, namely the Gauss linear scheme, due to its computational efficiency. Divergence terms are discretized using different schemes, where the linearUpwind scheme is applied for ∇∙U, the upwind scheme is applied for ∇∙k and ∇∙ω, and the Gauss linear scheme is applied for ∇∙ν_eff(∇U)^T. Laplacian terms are all discretized using the Gauss linear scheme with non-orthogonality correction. Interpolations are calculated using the linear scheme. Surface normal gradients are discretized using the non-orthogonality correction scheme, namely corrected. Wall distance is calculated using the meshWave method. It is noted that these spatial schemes are also used in the URANS simulations except that the time scheme is changed from steadyState to backward, which is a second-order implicit scheme for time discretization.

7. Solving Algorithms—fvSolution

The fvSolution file of RANS simulations in the example case is shown in Figure 8. The preconditioned conjugate gradient (PCG) linear solver [72] with the diagonal-based incomplete Cholesky (DIC) preconditioner is used to solve matrices of pressure and flux. Matrices of the rest parameters are solved by the preconditioned bi-conjugate gradient (PBiCGStab) linear solver [73] with the diagonal-based incomplete lower-upper (DILU) preconditioner. Both solvers have a solving tolerance of 1 × 10⁻¹¹. The SIMPLE algorithm [74] is used to perform the pressure-velocity coupling. In OpenFOAM, the simpleFoam solver can adopt both the SIMPLE algorithm and the SIMPLEC algorithm [75], where the former is selected by setting the consistent variable to a false value. Before running the simpleFoam solver, the potentialFoam solver is run for 10 iterations to initialize velocity, pressure, and flux fields. The gradient of the velocity field is cached to speed up the solver.

The fvSolution file of the URANS simulations in the example case is shown in Figure 9. For URANS simulations in the example case, the pimpleFoam solver is utilized, which adopts a hybrid algorithm of PISO [76] and SIMPLE. Configurations of linear solvers in the URANS simulation are similar to those in the RANS simulations, except that the relative tolerance of these linear solvers is set to be zero at the final iteration of the outer corrector for each parameter. The value of nCorrectors is set to be 2, which is within the recommended values (1 to 3) in the OpenFOAM user guide v9 [67]. Although the user guide also suggests that the nOuterCorrectors take a value from 50 to 1000, we use a larger value (2000) to ensure that the outer correction converges. The outer correction iterations work like the SIMPLE iterations, and so under-relaxation factors are applied.

8. Parallel Configurations—decomposeParDict

Simulations in this example case would likely have cell counts of over 30 million, and so running simulations in parallel is necessary. The SCOTCH algorithm [77] should be selected in the decomposeParDict file to perform mesh decomposition. For our research applications, simulations were performed on 128 CPU cores (Intel Xeon Gold 6148 processors running at 2.4 GHz) of the Irish Centre for High-End Computing (ICHEC) Kay server. For commercial applications, alternative processing capability would need to be sourced. The use of 128 CPU cores is not random but a decision based on a strong scalability test on the server where the RANS simulation is run for 6000 iterations with different numbers of CPU cores. Details of these tests are summarized in

Table 2. Results of the strong scalability test.

Number of CPU Cores	Wall-Clock Time (in h)	Speed-Up	Cells Per CPU Core
1	965.28 (estimated)	-	33,189,094
16	60.33	16	1,937,500
32	35.88	26.903	968,750
64	16.31	59.172	484,375
128	7.84	123.023	242,188
256	3.61	267.175	121,094

. Figure 8 shows the simulation time and the simulation speed-up. In an ideal condition, the simulation should scale linearly, i.e., using 64 CPU cores should speed up the simulation 64 times. However, Figure 10b shows that using 64 or fewer cores does not reach the idealized speed-up rate. On the contrary, a super-linear parallel speed-up is seen, where using 256 cores speeds up the simulation by a factor of 267 versus using 1 CPU core. This can be explained by the cache effects in parallel computing [78]. Similar super-linear parallel speed-up due to the cache effects was also reported in simulations by Cardiff et al. [79]. Although the computing time of the simulation using 256 cores is over 3 h shorter than using 128 cores, the queueing time on the server is significantly increased. Therefore, using 128 CPU cores provides a balance between the queueing time and the computing time for the example case.

9. Mesh Sensitivity Study

As mentioned in Section 3, the mesh quality can affect the accuracy of simulation results. Hence, a mesh sensitivity study should be performed. Typically, this involves creating a coarser mesh and a finer mesh, running the simulations for each, and comparing results. Ideally, the finer mesh will give results that are the same as the starting mesh. This indicates that there is no benefit to using the finer mesh and it would be wasteful of computational effort. Ideally, the coarser mesh will give different results, indicating that it is not accurate enough. In doing this, one can have confidence that the chosen mesh is appropriate for the task at hand. Figure 11 presents four meshes for our example case study to represent a mesh sensitivity study. The mesh described in Section 3 is represented here in Figure 11c.

The coarse mesh in Figure 11a has a background mesh where cells within have a length of 0.1 m. Based on this background mesh, refinement of six levels was applied, where two buffer layers were added between adjacent levels. The smallest cell within the coarse mesh has a thickness of 1.5625 × 10⁻³ m, which corresponds to a y⁺ value of approximately 10. The medium mesh (Figure 11b) has a background mesh that is one level finer than the coarse mesh. To better capture the flow characteristics in the region that is close to the geometry, a refinement box was used, where cells were firstly refined to level 4 and then refined to level 6 at the surface of the geometry. The smallest cell within the medium mesh has a thickness of approximately 7.8125 × 10⁻⁴ m, which corresponds to a y⁺ value of approximately 5. The finer mesh in Figure 11d has a similar scheme as the fine mesh (Figure 11c), except that it has a finer background mesh where cells within have a thickness of 0.03 m. Accordingly, the y⁺ value of the finer mesh (approximately 0.25) is smaller than that of the fine mesh (approximately 0.45).

Table 3. Cell counts of different meshes.

Configurations	Coarse Mesh	Medium Mesh	Fine Mesh	Finer Mesh
−10°	1,962,324	14,607,636	33,242,170	85,922,392
0°	1,973,227	14,575,020	33,189,094	86,008,513
10°	1,972,355	14,608,860	33,227,289	86,441,092

shows the cell counts of these meshes. It should be noted that with different configurations of the angle of attack, cell counts slightly vary for each mesh.

It is straightforward to perform a mesh sensitivity study in simulations using 2D structured meshes, where they can compare parameters calculated from simulations using meshes of different mesh density, i.e., different cell counts. However, for modelling wind effects on bridges where meshes should be 3D and unstructured, the cell count is no longer suitable to illustrate the mesh density. This is because meshes with a similar cell count can have very different cell sizes at the near-wall region, which is of interest in this type of study when using different refining schemes. Here, average cell sizes at the patch of interest (the bridge) should be used as an indicator of mesh density. Figure 12 demonstrates a mesh sensitivity study using the four meshes.

Aerodynamic coefficients are calculated using the presented four meshes and corresponding numerical configurations in the following sections. It can be found that the use of the coarse mesh and medium mesh can determine results that are up to over 80% different from the results determined using the finer mesh. In contrast, there is no appreciable change in results between using the fine mesh and the finer mesh where relative differences are below 2%. In addition, the use of the finer mesh requires more than three times the computational power than the fine mesh. Therefore, based on the results of the mesh sensitivity study, the fine mesh is deemed to be the suitable mesh for the study.

It should be noted that only four mesh schemes are presented here to provide readers with a sense of how a suitable mesh is developed using snappyHexMesh in OpenFOAM. However, in reality, it is likely that far more than four meshes would be developed and analyzed before choosing a final mesh.

10. Post Processing

10.1. Aerodynamic Force Coefficients

During each SIMPLE iteration and PISO or PIMPLE time step, aerodynamic forces, namely the drag force (F_x), lift force (F_z), and pitching moment (M_y), are calculated using the force function.

As this is an incompressible case, the air density (at 15 °C) is specified with a constant value of 1.225 kg/m³. The porosity variable is given a false value since porous effects are not considered in bridge deck aerodynamic studies. The CofR variable defines the center of rotation for calculating the pitching moment.

Output of the force function is stored in the postprocessing/directory, where both viscous and pressure parts of drag forces (F_x), lift forces (F_z), and pitching moments (M_y) are recorded at each time step and SIMPLE iteration in URANS simulations and RANS simulations, respectively. Summations of the viscous and pressure parts of each force should be taken to determine their full values. Then, aerodynamic coefficients of the bridge deck, namely the drag coefficient (C_Fx), lift coefficient (C_Fz), and moment coefficient (C_My), are calculated in Equations (5)–(7) based on aerodynamic forces:

C_Fx = 2 F_x/ρU² HL

(5)

C_Fz = 2 F_z/ρU² BL

(6)

C_My = 2 M_y/ρU² B^2 L

(7)

where F_x, F_z, M_y are drag force, lift force, and pitching moment, respectively; ρ is the density of the air, which is 1.225 kg/m³ and corresponds to the temperature at 15 °C; U is the reference velocity (P₁); B is the effective width of the bridge deck (P₄); H is the depth of the bridge deck (P₅); and L is the span of the bridge deck (P₆). Figure 13 shows the sample aerodynamic coefficients that would be calculated from the simulation results. These coefficients are essential to bridge design, especially in the topological optimization of bridge decks.

10.2. Visualization

Apart from the quantitative results that can be calculated using the various functions, simulation results can also be visualized. The authors recommended that this be done in ParaView [80]. There are many ways to visualize results. Figure 14a,b shows typical velocity plots that would be generated from our example case of RANS and URANS simulations, respectively. To illustrate the flow circulations, the line integral convolution (LIC) representation is used in Figure 15. The images clearly illustrate the simplifications of a steady-state analysis whereby the velocities appear smooth in Figure 14a and Figure 15a and more rotational in Figure 14b and Figure 15b. The URANS simulation better captures the intricacies of complex flow analysis. Figure 16 shows another parameter that bridge engineers are interested in, the wall shear stress on the bridge deck surface, where Figure 16a shows wall shear stress from a RANS simulation and Figure 16b shows that from a URANS simulation. Evidence of flow structures are present in Figure 16b that are not present in Figure 16a. Finally, using the Q-criterion, Figure 16a,b show the turbulence structures from RANS simulations and URANS simulations, respectively. The Q-criterion (Figure 17) defines vortices as areas where the vorticity magnitude is greater than the magnitude of the rate of strain.

These visualizations provide a powerful tool for bridge designers to help communicate results of wind effects on bridges. They offer an opportunity to visually demonstrate to architects, clients, planning authorities, and others involved in the design and commissioning process of new bridges the implications of one design over another. They offer a similarly useful tool for existing bridges where a variety of remediation works are being considered.

11. Validation of the Methodology

The methodology presented in this paper has been deliberately developed to be sufficiently generic so as to accommodate many types of long-span bridges while still being viable for industrial applications. It is designed to have enough detail for bridge designers to use that can also be trusted based on its implementation with two real bridges—the Queensferry Crossing in Edinburgh [22,23] and the Rose Fitzgerald Kennedy Bridge in Ireland [21,46].

Zhu et al. [23] simulated a full-scale three-dimensional CFD model of the middle span and central tower of the Queensferry Crossing that was created in OpenFOAM with the k-ω SST turbulence model. Atmospheric boundary layer inflows were configured based on wind profiles provided by a full-scale weather research and forecasting (WRF) model. CFD predictions were validated with field data, which were collected from an on-site structural health monitoring system. The simulated fluctuating wind field closely satisfied the characteristics of field data and demonstrated that the modelling approach had good potential to be used in practical bridge aerodynamic studies. Figure 18 below, reproduced from [23], shows wind velocities at a particular location estimated from CFD simulations with three different mesh refinements and verified with the field data recorded at the same location, indicating an average difference of just 10%.

Zhang et al. [46] developed a full-scale 3-dimensional CFD simulation model to replicate wind conditions at the Rose Fitzgerald Kennedy Bridge in Ireland. Wind velocities at four different points calculated by CFD simulations were compared with corresponding data collected from SHM field measurements. The calculated time-averaged wind velocities at four different locations on the bridge were shown to have relative differences of less than 10% to the measured wind velocities by anemometers 90% of the time. Figure 19, reproduced from [46], demonstrates these results for one of the sensors. In this example, the average relative difference within the comparison is 6%, and the underlying behavior within the time history of the wind velocity is generally predicted by the CFD simulations.

12. Conclusions

In this paper, a workflow is presented to develop CFD models in OpenFOAM to assess wind effects on long-span bridges that are computationally viable for industrial applications. Detailed configurations and corresponding considerations are described within each step. Comments on the validation of this methodology with field data collected from two-long span bridges are provided. By adopting this workflow, OpenFOAM can be used as an efficient tool that provides bridge engineers with key parameters in their design of long-span bridges against wind effects. In addition, the open-source nature makes this tool more affordable compared to commercial alternatives. Therefore, it can be foreseen that OpenFOAM will have an increasingly large community among civil engineering users.