Unsteady RANS Simulations of Flow around a Twin-Box Bridge Girder Cross Section

Abstract

The aerodynamic performance of bridge deck girders requires a thorough assessment and optimization in the design of long-span bridges. The present paper describes a numerical investigation of the aerodynamic characteristics of a twin-box bridge girder cross section in the range of angles of attack between −10.0° and +10.2°. The simulations are performed by solving 2D unsteady Reynolds-averaged Navier–Stokes (URANS) equations together with the k–ω shear stress transport (SST) turbulence model. The investigated Reynolds number (Re) based on the free stream velocity ($U_{\infty }$) and the height of the deck (D) is 31,000. The predicted aerodynamic characteristics such as the mean drag, lift and moment coefficients, are generally in good agreement with the results from the wind tunnel tests. Changes of flow patterns and aerodynamic forces with different angles of attack are investigated. Flow characteristics during one vortex shedding period are highlighted. Relative contributions of each of the two bridge decks to the overall drag and lift coefficients, with respect to the angle of attack, are also discussed.

1. Introduction

In wind-resistant bridge design, many factors should be considered, such as static wind load, wind forces due to turbulence, aerodynamic instability and vortex-shedding excitation. For bridge spans around 1 km and longer, twin-box girders are being increasingly used, due to their favorable aerodynamic properties in term of flutter stability. This is due to a favorable effect of the central gap on the surface pressure distribution, and a higher ratio between the twisting and the heaving eigen-frequencies [1,2]. On the other hand, this configuration is prone to vortex-induced vibration, since vortices shed from the windward box can cause a significant excitation of the leeward box and thereby important oscillations of the bridge girder [3,4]. In the early design stage of long-span bridges, several bridge deck design alternatives are usually considered, and their aerodynamic performance must be assessed to ensure the most feasible design for the specific project.

Both experimental and numerical approaches are available for the assessment of the bridge deck aerodynamic performance [4,5]. Extensive wind tunnel experiments are normally performed to optimize the design of long-span bridges. However, wind tunnel investigations have limitations in terms of time, cost, applicable Reynolds numbers and the representation of the details such as railings and stabilizers at small scales. For an improved bridge design, computational fluid dynamics (CFD) based simulations can facilitate a deeper understanding of complex flow conditions around a bridge girder, as demonstrated by [6,7,8,9,10,11,12,13].

The application of a 2D unsteady Reynolds-averaged Navier–Stokes (URANS) model for flow around simple rectangular geometry has been validated thoroughly in previous research [14,15,16,17,18]. Mannini et al. [16] applied 2D URANS equations with advanced turbulence modeling closures for predicting flow around a 5:1 rectangular cylinder at Re = 1 × 10⁵. They found that the predicted force coefficients are in reasonable agreement with the experimental data reported by Schewe [19,20]. Ong [17] conducted CFD simulation for the same structure at high Re numbers (5 × 10⁵, 1 × 10⁶, 1.5 × 10⁶ and 2 × 10⁶) using the 2D URANS together with the standard high Reynolds number k–ω model. For engineering design purposes, it gave reasonably good agreements with the published experimental data in terms of time-averaged drag coefficients, Strouhal numbers and time-averaged base pressure coefficient. Patruno et al. [18] studied the flow field around the 5:1 rectangular cylinder with an angle of attack of 0 to 4 degrees using both Large Eddy Simulation (LES) and URANS (i.e., k–ω shear stress transport (SST)) turbulence models. They found that LES and URANS approaches appeared to provide comparable results in terms of accuracy of first and second order pressure statistics measured at the central section of the considered prism. The asymmetric characteristics of the flow field induced by small angles of attack were under-estimated by LES as compared to URANS simulations, and the experimental data were between them generally.

The main goal of the present study is to evaluate the validity of 2D URANS simulations with k–ω SST turbulence model [21] for the flow around a twin-box bridge girder cross-section. One of the preliminary bare deck designs of the Halsafjord suspension bridge in Norway will be considered in the present study, see Figure 1. This bare deck design was chosen due to the availability of experimental data for validation studies. Several CFD studies of flow characteristics around the sharp-edged twin-boxes, particularly with asymmetric geometry as in the present case, have been carried out [22,23,24]. These studies have mainly dealt with the overall flow conditions and their variations with the relative gap size, Reynolds number and turbulence models used in the simulations. To the authors’ knowledge, there is a limited number of studies addressing the flow characteristics around the bridge girder subject to flow at different angles of attack. The present study will provide a thorough 2D CFD investigation of the flow conditions around the twin-deck inclined to the approaching flow at angles of attack (AoA) ranging from −10.0 to +10.2 degrees (°).

The study utilizes data from a wind tunnel investigation for validation. The tests were performed with a 1:40 section model at Re = 31,000, as reported in [25]. The cross-section geometry and the Reynolds number in the simulations are thus adopted as for the section model, which was exposed to a uniform flow. The results are discussed in terms of the time-averaged drag coefficient ($\overline{C_D}$), the time-averaged lift coefficient ($\overline{C_L}$) and the time-averaged moment coefficient ($\overline{C_M}$) in the range of AoA from −10.0° to +10.2°, respectively. The time histories of the force coefficients as well as the flow physics changes with different AoA are investigated. The vortex shedding process and the associated variation of the instantaneous drag and lift coefficients for each of the two decks with different AoA are also presented and discussed in detail.

2. Mathematical Formulation

2.1. Flow Model

The Reynolds-averaged equations for conservation of mass and momentum are given by:

$$\frac{\partial u_j}{\partial x_j}=0$$

(1)

$$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+v\frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}-\frac{\partial \overline{u_i^{\text{'}}{u}_j^{\text{'}}} }{\partial x_j}$$

(2)

where i, j = 1, 2. Here $x_1$ and $x_2$ denote the streamwise and cross-stream directions respectively; $u_1$ and $u_2$ are the corresponding mean velocity components; $\overline{u_i^{\text{'}}{u}_j^{\text{'}}}$ is the Reynolds stress component, where $u_i^{\text{'}}$ denotes the fluctuating part of the velocity; $p$ is the pressure; $\rho$ is the density of the fluid and $t$ is the time.

The k–$\omega$ SST turbulence model [26] used in the present study is a blending of the k–$\omega$ and the k–$\epsilon$ models. The original k–$\omega$ model of [27] is implemented in the near-wall region and the standard k–$\epsilon$ model of [28] in the outer wake region and in the free shear layers. Following [26], the equations for the k–$\omega$ SST turbulence model is taken as:

$$\frac{D\left(\rho k\right)}{Dt}=\widetilde{P_k}-{\beta }^{*}\rho \omega k+\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_k{\mu }_t\right)\frac{\partial k}{\partial x_j})\right]$$

(3)

$$\frac{D\left(\rho \omega \right)}{Dt}=\alpha \rho S^2-\beta \rho {\omega }^2+\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\frac{\partial \omega }{\partial x_j})\right]+2\left(1-F_1\right)\rho {\sigma }_{\omega 2} \frac{1}{\omega } \frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}$$

(4)

where $\widetilde{P_k}$ is given by:

$$\widetilde{P_k}=\min \left[{\mu }_t\frac{\partial u_i}{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right), 10{\beta }^{*}\rho k\omega \right]$$

(5)

Here ${\varphi }_1$ represents any constant in the original k–$\omega$ turbulence model (i.e., ${\sigma }_{k1},{\beta }_1,{\sigma }_{\omega 1})$ and ${\varphi }_2$ represents any constant in the original $k–\epsilon$ turbulence model (i.e., ${\sigma }_{k2},{\beta }_2,{\sigma }_{\omega 2})$.$\varphi$, the corresponding constant of the new $k-\omega$ SST turbulence model given by Equations (3) and (4), is

$$\varphi =F_1{\varphi }_1+\left(1-F_1\right){\varphi }_2$$

(6)

$$F_1=\tan \mathrm{h}({\arg }_1{}^4)$$

(7)

$${\arg }_1=min\left[\max \left(\frac{\sqrt{k}}{0.09\omega y},\frac{500\nu }{y^2\omega }\right),\frac{4\rho {\sigma }_{\omega 2}k}{C{D}_{k\omega }{y}^2}\right]$$

(8)

$$C{D}_{k\omega }=max\left(2\rho {\sigma }_{\omega 2}\frac{1}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j},{10}^{-10}\right)$$

(9)

Here $y$ is the distance to the nearest wall and $C{D}_{k\omega }$ is the positive portion of the cross-diffusion term of Equation (4).

The turbulence viscosity can be estimated by:

$${\nu }_t=\frac{a_{1}k}{\max \left(a_1\omega ,\mathrm{S}{F}_2\right)}$$

(10)

where S is the invariant measure of the strain rate and $F_2$ is given by:

$$F_2=\tan \mathrm{h}({\arg }_2{}^2),{\arg }_2=\max \left(\frac{2\sqrt{k}}{0.09\omega y},\frac{500\nu }{y^2\omega }\right)$$

(11)

The empirical constants of the k–$\omega$ SST model are [19]:

${\beta }^{*}=0.09,a_1=0.31,{\alpha }_1=0.5532,{\alpha }_2=0.4403,{\beta }_1=0.075,{\beta }_2=0.0828,{\sigma }_{k1}=0.85034,{\sigma }_{k2}=1.0$ and ${\sigma }_{\omega 1}=0.5,{\sigma }_{\omega 2}=0.856$16.

2.2. Numerical Simulation Scheme, Computational Domain and Boundary Conditions

The open source CFD code OpenFOAM [29] is used in the present study. The pressure implicit with splitting of operators (PISO) scheme is chosen for the solver. The spatial schemes for gradient, Laplacian and divergence are Gauss linear, Gauss linear corrected and Gauss linear schemes, respectively. All these schemes are in second order accuracy. The second order Crank–Nicolson scheme is used for the time integration.

Figure 2 shows the investigated twin-box bridge girder cross section and the corresponding computational domain. The size of the entire computational domain is 6 m by 1.55 m, which is equivalent to 91.0D by 24.8D, where D is the height of the girder cross section. The origin of the coordinate system is located at the left lower corner of the domain. The flow inlet boundary is 18.8D upstream from the center of the gap between the decks and the flow outlet boundary is 77.2D downstream from the center of the gap. The top and bottom boundaries are positioned at a distance of 12.4D from the center of the gap between the decks. Although the shear center of the bridge deck section does not coincide with the mid height of the deck, the rotation center is adopted in the middle of the deck height. This facilitates a direct comparison with the experimental results which also refer to such a rotational center. The deck position and the domain size in the Y direction are same as in the wind tunnel working section. The boundary conditions used for the present numerical simulations are as follows:

• A uniform flow, $u_1=U_{\infty }$, $u_2=0$, is set at the inlet boundary; the pressure is specified as zero normal gradient at the inlet boundary. $k$ and $\omega$ at the inlet boundary are set equal to:

  $$k_{inlet}=0.5{\left(U_{\infty }{I}_u+U_{\infty }{I}_v+U_{\infty }{I}_w\right)}^2$$

  (12)

  where $I_u$ and $I_v$ are the turbulence intensities in X and Y directions, respectively. $I_u=12\%$ and $I_v=8\%$ are taken from [25].

  $${\omega }_{inlet}=k_{inlet}{}^{0.5}/\left(C_{\mu }{}^{0.25}l\right)$$

  (13)

  where $C_{\mu }=0.09$ is the empirical constant specified in the turbulent model and turbulence length $l$ = 0.07D [30,31,32]. Effects of $l$ on the calculated results have been studied by using a much lower value $l$ = 0.04D, and small variations (less than 0.27%) are observed in the aerodynamic quantities, i.e., mean drag coefficient, root mean square (r.m.s.) of the lift coefficient and Strouhal number.
• Along the outlet boundary, $u_1, u_2$, $k$ and $\omega$ are specified as zero normal gradient. The pressure is specified as zero. The zero pressure outlet boundary condition has been widely used to calculate the unsteady flow around bluff bodies [33,34,35]. The distance used in the present study for the downstream 77.2D is considerably longer than the value of 50D previously used by [10]. It is considered that the effect of the outlet boundary condition on the numerical results is negligible.
• On the deck surfaces, no-slip boundary condition is specified (i.e., $u_1=u_2=0$). The pressure is set as zero normal gradient. k is fixed at 0 and $\omega$ is calculated as follows [26]:

  $${\omega }_{deck}=10\times \frac{6\nu }{{\beta }_1{\left(\Delta y_1\right)}^2}$$

  (14)
• For the top and bottom boundaries, $u_1=u_2=0,$ the pressure is set as zero normal gradient; k is fixed at 0. $\omega$ is specified as zero normal gradient.

2.3. Grid and Time Resolution Tests

Convergence studies considering different grid sizes and the time steps have been carried out for AoA ranging from −10.0° to +10.2°. Time-averaged drag coefficient ($\overline{C_D}$), time-averaged lift coefficient ($\overline{C_L}$), time-averaged moment coefficient ($\overline{C_M}$) as well as Strouhal number (St) are considered in the tests. The time-averaged values are obtained with a duration of 30 vortex shedding cycles after the numerical results have repeated their cycles. Here the drag coefficient ($C_D$), lift coefficient ($C_L$), moment coefficient ($C_M$) and Strouhal number (St) are defined as follows:

$$C_D=\frac{F_{drag}}{\frac{1}{2}\rho U_{\infty }{}^{2}D}$$

(15)

$$C_L=\frac{F_{lift}}{\frac{1}{2}\rho U_{\infty }{}^{2}b}$$

(16)

$$C_M=\frac{M}{\frac{1}{2}\rho U_{\infty }{}^2{b}^2 }$$

(17)

$$St=\frac{fD}{U_{\infty } }$$

(18)

where $F_{drag}$, $F_{lift}$ are the along-wind and the cross-wind force components acting on the decks per unit length and b is the width of the two decks, i.e., of the solid part of the cross-section. $M$ is the overturning moment about the aforementioned rotational center per unit length, and is positive in the clockwise direction. The frequency of vortex shedding f is obtained from C_L. Figure 3 shows the sign convention. Like the overturning moment, AoA is defined positive in the clock-wise (nose-up) direction in the cases of wind coming from left to right.

The results of both the grid size and the time step convergence studies in terms of the force coefficients are shown in

Table 1. Aerodynamic properties for different grid and time resolutions for AoA > 0°.

$\mathit{A}\mathit{o}\mathit{A}.\left(\text{°}\right)$	Case	Elements	Δt (s)	$\overline{{\mathit{C}}_{\mathit{D}}}$	$\overline{{\mathit{C}}_{\mathit{L}}}$	$\overline{{\mathit{C}}_{\mathit{M}}}$	St	$\frac{\mathbf{\Delta }}{\overline{{\mathit{C}}_{\mathit{D}}}}(\%)$	$\frac{\mathbf{\Delta }}{\overline{{\mathit{C}}_{\mathit{L}}}}(\%)$	$\frac{\mathbf{\Delta }}{\overline{{\mathit{C}}_{\mathit{M}}}}(\%)$	$\frac{\mathbf{\Delta }}{\mathit{S}\mathit{t}}(\%)$
0.1	M1	138,354	5.00 × 10−5	1.164	−0.154	0.077	0.224	-	-	-	-
0.1	M2	198,834	5.00 × 10−5	1.129	−0.176	0.078	0.216	−3.0	−14.0	1.3	−3.7
0.1	M3 *	288,034	5.00 × 10−5	1.125	−0.178	0.075	0.216	−0.3	−1.7	−2.8	0.0
0.1	M3T1	288,034	2.50 × 10−5	1.142	−0.180	0.082	0.225	1.5	−0.6	8.2	4.3
1.5	M1	138,354	5.00 × 10−5	1.088	−0.051	0.127	0.203	-	-	-	-
1.5	M2	198,834	5.00 × 10−5	1.096	−0.057	0.125	0.201	0.7	−12.3	−1.8	−0.7
1.5	M3 *	288,034	5.00 × 10−5	1.099	−0.057	0.123	0.205	0.2	−0.7	−1.8	2.0
1.5	M3T1	288,034	2.50 × 10−5	1.114	−0.047	0.126	0.196	1.4	18.2	2.2	−4.4
3.2	M1	138,354	5.00 × 10−5	1.080	0.141	0.182	0.201	-	-	-	-
3.2	M2	198,834	5.00 × 10−5	1.087	0.160	0.181	0.191	0.7	14.0	−0.5	−5.1
3.2	M3 *	288,034	5.00 × 10−5	1.091	0.154	0.180	0.192	0.3	−3.8	−0.8	0.4
3.2	M3T1	288,034	2.50 × 10−5	1.101	0.147	0.180	0.200	1.0	−5.0	−0.2	4.4
4.4	M1	138,354	5.00 × 10−5	1.106	0.245	0.218	0.186	-	-	-	-
4.4	M2 *	288,034	5.00 × 10−5	1.106	0.228	0.210	0.192	0.0	−6.8	−3.8	3.4
4.4	M3	545,954	5.00 × 10−5	1.083	0.238	0.211	0.200	−2.1	4.1	0.7	4.0
4.4	M2T1	545,954	2.50 × 10−5	1.112	0.241	0.211	0.196	2.7	1.5	−0.1	−2.2
6.3	M1	138,354	2.50 × 10−5	1.429	0.260	0.177	0.261	-	-	-	-
6.3	M2 *	288,034	2.50 × 10−5	1.111	0.313	0.235	0.250	−22.3	20.4	32.7	−4.1
6.3	M3	545,954	2.50 × 10−5	1.096	0.313	0.239	0.240	−1.4	0.0	1.5	−3.9
6.3	M2T1	545,954	1.25 × 10−	1.107	0.325	0.240	0.244	1.0	3.7	0.5	1.6
8.1	M1	138,354	2.50 × 10−5	1.820	0.338	0.207	0.126	-	-	-	-
8.1	M2	198,834	2.50 × 10−5	1.654	0.371	0.202	0.136	−9.1	9.8	−2.4	8.5
8.1	M3 *	288,034	2.50 × 10−5	1.628	0.387	0.206	0.142	−1.6	4.4	2.0	4.2
8.1	M3T1	288,034	1.25 × 10−	1.636	0.395	0.207	0.137	0.5	2.0	0.7	−3.3
10.2	M1	198,834	2.50 × 10−5	2.045	0.263	0.148	0.125	-	-	-	-
10.2	M2	272,914	2.50 × 10−5	2.037	0.210	0.141	0.107	−0.4	−20.1	−5.0	−14.6
10.2	M3 *	362,514	2.50 × 10−5	2.039	0.208	0.132	0.110	0.1	−1.0	−5.8	2.8
10.2	M3T1	362,514	1.25 × 10−5	2.057	0.200	0.130	0.101	0.9	−4.0	−1.8	−7.9

In the columns 9–12, ∆ indicates the variation of the aerodynamic parameter relative to the previous case.

and

Table 2. Aerodynamic parameters for different grid and time resolutions for AoA < 0°.

AoA (°)	Case	Elements	Δt (s)	$\overline{{\mathit{C}}_{\mathit{D}}}$	$\overline{{\mathit{C}}_{\mathit{L}}}$	$\overline{{\mathit{C}}_{\mathit{M}}}$	St	$\frac{\mathbf{\Delta }}{\overline{{\mathit{C}}_{\mathit{D}}}}(\%)$	$\frac{\mathbf{\Delta }}{\overline{{\mathit{C}}_{\mathit{L}}}}(\%)$	$\frac{\mathbf{\Delta }}{\overline{{\mathit{C}}_{\mathit{M}}}}(\%)$	$\frac{\mathbf{\Delta }}{\mathit{S}\mathit{t}}(\%)$
−1.4	M1	138,354	5.00 × 10−5	1.277	−0.190	0.024	0.228	-	-	-	-
−1.4	M2	198,834	5.00 × 10−5	1.230	−0.208	0.020	0.214	−3.7	−9.2	13.6	−6.3
−1.4	M3 *	288,034	5.00 × 10−5	1.230	−0.209	0.019	0.217	0.0	−0.7	4.9	1.5
−1.4	M3T1	288,034	2.50 × 10−5	1.205	−0.210	0.024	0.211	−2.0	−0.3	−22.2	−2.7
−2.9	M1	138,354	5.00 × 10−5	1.345	−0.208	−0.053	0.226	-	-	-	-
−2.9	M2	198,834	5.00 × 10−5	1.339	−0.207	−0.038	0.218	−0.5	0.0	27.2	−3.8
−2.9	M3 *	288,034	5.00 × 10−5	1.340	−0.200	−0.038	0.220	0.1	3.7	1.0	1.0
−2.9	M3T1	288,034	2.50 × 10−5	1.318	−0.191	−0.039	0.226	−1.6	4.6	−2.6	2.6
−4.0	M1	138,354	5.00 × 10−5	1.512	−0.251	−0.080	0.212	-	-	-	-
−4.0	M2	198,834	5.00 × 10−5	1.482	−0.219	−0.066	0.215	−2.0	12.9	17.2	1.3
−4.0	M3 *	288,034	5.00 × 10−5	1.487	−0.217	−0.066	0.212	0.4	0.9	0.6	−1.5
−4.0	M3T1	288,034	2.50 × 10−5	1.485	−0.216	−0.065	0.213	−0.1	0.2	0.9	0.4
−6.2	M1	138,354	5.00 × 10−5	1.512	−0.354	−0.165	0.203	-	-	-	-
−6.2	M2	198,834	5.00 × 10−5	1.517	−0.362	−0.163	0.196	0.3	−2.3	−1.5	−3.6
−6.2	M3 *	288,034	5.00 × 10−5	1.519	−0.364	−0.163	0.198	0.1	−0.6	0.1	1.2
−6.2	M3T1	288,034	2.50 × 10−5	1.512	−0.381	−0.161	0.203	−0.5	−4.6	−0.8	2.1
−8.1	M1	138,354	1.25 × 10−5	2.094	−0.681	−0.181	0.160	-	-	-	-
−8.1	M2	198,834	1.25 × 10−5	1.963	−0.653	−0.188	0.180	−6.3	4.2	3.5	12.8
−8.1	M3 *	288,034	1.25 × 10−5	2.023	−0.664	−0.182	0.182	3.1	−1.8	−3.3	1.1
−8.1	M3T1	288,034	6.25 × 10−6	2.064	−0.649	−0.191	0.176	2.0	2.2	5.3	−3.2
−10.0	M1	138,354	2.50 × 10−5	2.567	−0.793	−0.135	0.182	-	-	-	-
−10.0	M2	198,834	2.50 × 10−5	2.473	−0.733	−0.185	0.100	3.6	7.6	37.1	−45.1
−10.0	M3 *	288,034	2.50 × 10−5	2.526	−0.733	−0.176	0.104	−2.2	−0.1	−5.1	4.0
−10.0	M3T1	288,034	1.25 × 10−5	2.507	−0.727	−0.173	0.095	0.7	0.9	−1.8	−8.9

In the columns 9–12, ∆ indicates the variation of the aerodynamic parameter relative to the previous case.

, where $\Delta t$. represents the time step in seconds. The relative change of results between two consecutive meshes or time steps is investigated in terms of $\overline{ C_D}$, $\overline{C_L}$, $\overline{C_M}$ as well as St. The cases with the converged solution in terms of grid size and time-step are marked with ‘*’ and, adopted for the further investigations.

For most of the cases, the grid resolution tests have been carried out with 3 sets of meshes with 138,354 (M1), 198,834 (M2) and 288,804 (M3) elements and two different time-steps (i.e., M3 with time-step of 5.00 × 10⁻⁵ s, and M3T1 with time-step of 2.50 × 10⁻⁵ s), see Table 1 and Table 2. For example, in the AoA = +0.1° case, the relative change of the $\overline{C_D}$ values between M1 and M2 meshes is 14.0%. The difference decreases to 1.7% when the mesh is further refined from M2 to M3, i.e., a negligible change is observed. A similar trend is observed for the other cases. For the cases of AoA = +4.4°, +6.3° and +10.2°, 138,354 to 545,954 (or 362,514) elements were used for the grid resolution tests. The convergence is considered to be achieved when the relative variation of the aerodynamic quantities between two cases is lower than 5% or if the parameters change by less than 0.01 in magnitude. For instance, in the AoA = +0.1° case, the relative change of $\overline{C_M}$ values between M3 and M3T1 cases is 8.2%, the absolute change of the $\overline{C_M}$ values is less than 0.01, see Table 1. An example of the mesh structure with 288,804 elements, in the case of AoA = +0.1°, is presented in Figure 4.

For the AoA = +0.1° case, the values of $y^{+}$ over the entire wall boundaries of two decks range from 0 to 7.8 with an average value from 1.0 to 1.5, where $y^{+}=u_{*}\Delta y_1/\nu =$ $\sqrt{{\tau }_w/\rho }·\Delta y_1/\nu$, where $u_{*}$ and ${\tau }_w$ denote the wall friction velocity and wall shear stress, respectively. The maximum $y^{+}$ value is found at the windward deck at the leading edge lower corner, where high flow acceleration exists. A similar procedure is taken to calculate the averaged $\overline{y^{+}}$ for other AoA cases, and the averaged values of $\overline{y^{+}}$ give a range of 1.0 to 2.1. All refinement studies of the flow at different AoA are carried out with the fixed $y^{+}$ values in order to minimize influence of the boundary layer calculation.

Overall, it is concluded that the present simulations marked with ‘*’ in Table 1 and Table 2 can provide satisfactory spatial and time resolutions for all the different AoA cases at Re $=$ 31,000.

3. Results and Discussion

3.1. Validation of Numerical Model

To assess the validity of URANS with the k–$\omega$ SST turbulence model for flow around twin-box girder bridge cross section, a comparison has been carried out between the present simulations and the wind tunnel test results from [18]. The time-averaged aerodynamic quantities such as $\overline{C_D}$, $\overline{C_L}$ and $\overline{C_M}$, in the range of AoA from −10.0° to +10.2° at Re $=$ 31,000 are studied, see Figure 5. A good overall agreement between the present numerical simulations results and the experimental data is observed. Generally, the presently predicted $\overline{C_D}$ are in good agreement with the experimental data for all AoA. However, a relatively larger deviation from the experimental data is found at AoA = +6.3° and +8.1°. This may be due to a change in flow physics, which will be discussed in detail in Section 3.3.

The numerically based $\overline{C_L}$ and $\overline{C_M}$ values are generally in good agreement with the experimental counterparts. However, significant reductions of $\overline{C_L}$ and $\overline{C_M}$ at AoA = +10.2° are observed, where $\overline{C_L}$ is under-predicted by 57% and $\overline{C_M}$ by 31% as compared to the experimental data. Such a reduction in magnitude of $\overline{C_L}$ and $\overline{C_M}$ at large AoA are also reported in Nieto et al. [8], who performed the 2D simulations of the Stonecutters bridge deck cross section with the k–$\omega$ SST turbulence model. A similar over- or under-prediction of $\overline{C_L}$ and $\overline{C_M}$ at high AoA by the 2D URANS model can be found in several previous studies [3,6,7,8,9,13]. This might be attributed to the fact that the present 2D numerical model cannot account for the fluctuating flow in the axial direction (Z-direction), while the 3D flow effects become significant with the increase of AoA. The higher energy dissipation in 3D cannot be correctly evaluated in 2D, i.e., larger pressure fluctuations and stronger vortex generation on the body surface are observed in 2D analysis when compared to 3D analysis, see also [36]. In the normal operation condition, the bridge is mainly exposed to a relatively small AoA ($<3°)$, associated with the bridge motion and large-scale turbulence. For this range of angles, the 2D simulation can provide good and quick results, and can thus be an efficient evaluation tool for the bridge design, particularly as a screening tool for the bridge girder designs.

The predicted St, shown in Figure 5d, does not vary significantly (from 0.18 to 0.22) for AoA ranging from −8.1° to +4.4°. For the case of AoA = +6.3°, St increases to 0.250. This may be due to the flow structure change associated with an interaction of the vortices on the different parts of the deck surface, which will be discussed in detail in Section 3.3. For the cases of AoA = +8.1° and +10.2°, St gradually decreases to 0.142 and 0.110, respectively. The decrease in St at large AoA is consistent with an increase of the projected area in the cross-flow direction.

Overall, it appears that the present 2D RANS simulations with the k–$\omega$ SST turbulence model are in satisfactory agreement with the experimental data under the same flow conditions, especially at small AoA. The discrepancies of the aerodynamic coefficients between numerical and experimental results at large AoA (i.e., +10.2°) may be due to the three-dimensional effects of the flow structure along the spanwise direction.

3.2. Vortex Formation around the Decks in One Vortex Shedding Period

One of the main concerns in bridge design is the vortex-induced vibration (VIV) behavior of bridge girders. In order to investigate the vortex formation around the decks, eight time instants during one cycle of the C_L fluctuation at AoA = +0.1° are studied, as defined in Figure 6. It is observed that the values of the amplitudes of the maximum and minimum $C_L$ are not the same, and this is mainly due to the asymmetric geometry of the decks.

Figure 7 shows the vorticity contour plots for the eight time instants for AoA = +0.1°. During one C_L fluctuating period, it is observed that the vortices shed from the downstream upper and lower corners of the windward deck. Then the vortices propagate through the gap between the two decks and attack the leeward deck periodically. At the time instant 1, where C_L is the maximum, the lower vortex dominates the flow at the downstream edge of the windward deck what generates the lift force. Until the time instant 5, the upper vortex at the windward deck increases in size, and reduces the lift force gradually. As the vortices move to the leeward deck, the new vortices are generated from the downstream vertical edge of the windward deck. The upper and the lower vortices are different in size due to the asymmetric geometry of the decks. Vortices shed from the leeward edge of the downstream deck are very limited in size and strength due to a more streamlined form of this part of the cross-section.

3.3. Flow Characteristics at Different Angles of Attack (AoA)

Figure 8 shows the time histories of C_D, C_L and C_M of the AoA = +0.1° case, respectively. The normalized forces on each deck, as well as the sum of the normalized forces on the two decks are shown. It is observed that C_L and C_M vary at a more or less single frequency. The total drag coefficient C_D also includes a frequency component at twice the frequency of the vortex shedding. The values of C_D from the individual boxes do not follow a pure harmonic function. This is mainly attributed to the generation and the downstream action of the asymmetric vortices, see Figure 7. Figure 9 shows the time history of C_D for different AoA, i.e., +3.2°, +6.3°, −6.2° and −8.1°. By comparing Figure 8a and Figure 9a, it can be observed that the time history of C_D for AoA = +0.1° has different trend from that for AoA = +3.2°. This is due to the different vortex formation between these two cases. However, the time histories of C_L and C_M show similar behavior in both cases. By comparing Figure 9a–d, it is found that the time history of C_D shows significant changes for different AoA. The phase differences between the force coefficients on the two decks can be explained by the vortex shedding variations in Figure 7, i.e., the vortex shedding is generated first at the windward deck, and is reattached at the leeward deck.

Figure 9b shows time history of C_D for the AoA = +6.3° case. At this AoA, the flow characteristics begin to change. The second peak of C_D in one period becomes smaller as AoA increases. At AoA = +6.3°, the time history of C_D becomes a single harmonic, see Figure 9b. As the AoA increases, the projected area in the cross-flow direction increases and the drag force fluctuation of the windward deck decreases. The standard deviation (RSD) of C_D relative to its mean value, is less than 0.6% at the windward deck at AoA = +6.3°. It means that the drag force at the windward deck is almost stable, see Figure 9b. To study flow pattern change in detail, the C_L time histories for different AoA are looked into, as shown in Figure 10. For the AoA = +6.3° case, there are more vortex cycles than the other cases in the same time interval, corresponding to a higher St number value at this angle, see Figure 10b. Unlike the other cases, C_L at the windward deck is stable with little fluctuations during period at AoA = +6.3°. Figure 11a,b show instantaneous streamlines of the AoA = +6.3° case at the time instants of the maximum C_L at $t{U}_{\infty }/D$ = 304 and the minimum C_L at $t{U}_{\infty }/D$ = 368, respectively. It appears that the flow structure at the windward deck does not change significantly at these two time instants by comparing Figure 11a,b. However, the flow patterns at the leeward deck change significantly between these time instants. In particular, the vortices at the upper surface of the leeward deck and at the trailing edge of the leeward deck are observed to have strong interaction shown in Figure 11a,b. At the leeward deck, the interaction of the vortices on the upper deck surfaces seems to result in a flow separation over the entire surface.

Figure 12 show instantaneous streamlines of the AoA = +10.2°case at $t{U}_{\infty }/D$ = 504. At AoA = +10.2°, the vortices at the upper part of both the windward and the leeward decks encompass the entire surfaces. The flow separation leads to the loss of lift as indicated by a reduced C_L value, see Figure 5b.

In the range of negative AoA from −1.5° to −6.2°, the time histories of C_D, C_L and C_M have similar characteristics as those in the range of AoA from +0.1° to +4.4°. Figure 9c shows the time history of C_D for AoA = −6.2°, where a double peak within the main C_D cycle can be observed. Figure 9d shows the time history of C_D for the case of AoA = −8.1°. For the AoA = −8.1° case, the variation in the force coefficients is more random than that for the AoA = −6.2° case, see Figure 9 and Figure 10. Unlike the other cases depicted above, the time histories of aerodynamic quantities are no longer periodic. Generally, the time history of C_D shows a more random behavior than C_L by comparing Figure 9 and Figure 10. For the cases of AoA = −6.2° and AoA = +3.2°, the total C_D has a double peak due to the vortices generation at the windward deck, see Figure 9a,c. For the AoA = +6.3° case, the total C_D has one peak per period because the flow structure at the windward deck in one vortex shedding cycle changes significantly less than at the leeward deck. This explains the small variation of C_D at the windward deck observed in Figure 9b.

Figure 13a–c show instantaneous streamlines for the AoA = −10.0° case. For the AoA = −10.0° case, the flow both at the windward deck and the leeward deck changes significantly compared to the AoA = +0.1° case (small angle) shown in Figure 7. Thus, the different time-histories of the force coefficients reflect the differences in the flow pattern. For the AoA = +10.2° case in Figure 12, the big vortex at the upper surface of the leeward deck is nearly stationary and only small changes of the vortex structure are observed during the period. However, for the AoA = −10.0° case shown in Figure 13a–c, the vortices generated at lower surface of the leeward deck are not stationary. The vortex formation and the flow pattern change significantly. The lower inclined surface of the windward deck for the AoA = −10.0° case has a smaller angle to the flow. Thus, a weaker vortex generation mechanism is observed as compared to the conditions at the upper deck surface of the AoA = +10.2° case shown in Figure 12.

3.4. Contribution of Each Deck to $\overline{C_D}$ and $\overline{C_L}$

The relative contribution of the windward and the leeward deck, to the overall $\overline{C_D}$ and $\overline{C_L}$ is further presented in Figure 14. The $\overline{C_L}$ contributions are calculated based on the absolute mean values of each deck.

At AoA = +0.1°, the contribution of the windward deck to C_D is over 80%, and is less than 20% for the leeward deck, i.e., the along-wind force is mainly applied to the windward deck. As AoA increases, the contribution of the leeward deck to $\overline{C_D}$ increases. This is because the projected area of the leeward deck in the cross-flow direction increases as AoA increases, and the shielding effect of the upstream deck is reduced. The contribution of each deck to $\overline{C_D}$ is different for all the cases. For instance, the contributions of the two decks to $\overline{C_D}$ at AoA = +10.2° are nearly equal, i.e., 50% from the windward deck and 50% from the leeward deck. However, for the case of AoA = −10.0°, 30% contribution to $\overline{C_D}$ is from the windward deck, and 70% is from the leeward deck. The cross-sections of the bridge girders have asymmetric geometry, which result in different flow structures in the cases with same magnitude of AoA but in opposite signs.

Figure 12 shows the streamlines of the case of AoA = +10.2°, where the accelerated flow at the upper surface of the windward deck does not attack the leeward deck. Figure 13a–c show the streamlines of the AoA = −10.0° case, the accelerated flow at the upper surface of the windward directly hits the vertical edge of the leeward deck. Thus, even though two cases have same magnitude of AoA, the flow pattern and the contributions of each deck to the aerodynamic quantities are different. At all the negative AoA the leeward deck contributes more $\overline{C_D}$ than at the corresponding positive AoA. This is due to the accelerated flow at the windward deck acting on the leeward deck. For the positive AoA cases, the accelerated flow at the windward deck passes by the leeward deck and does not hit the leeward deck. The vortex at the upper surface of the windward deck prevents the interaction, see Figure 12.

For the contribution of each deck to $\overline{C_L}$, it is observed that the contribution of the windward deck to $\overline{C_L}$ increases as AoA increases. At AoA = +8.1°, nearly 100% lift force is generated by the windward deck and almost no contribution from the leeward deck. Then, at AoA = +10.2° the $\overline{C_L}$ contribution from the windward deck drops to around 90%, where the abrupt lift drop is observed. A similar trend can be observed in the negative AoA. From AoA = −1.4° to −6.2°, the contribution of the windward deck to $\overline{C_L}$ increases. For the cases of AoA = −8.1° and −10.0°, the contribution of the windward deck to $\overline{C_L}$ drops. The latter two cases have shown that the flow pattern begins to change and becomes chaotic. Thus, this C_D and C_L contribution study can offer a general overview of the aerodynamic responses and flow physics of the bridge. In particular, it is a good qualitative tool to figure out where the flow pattern begins to change.

Figure 15 and Figure 16 show the variations of the relative standard deviation (RSD) of C_D and C_L to each deck versus AoA, respectively. The RSD of C_D and C_L are obtained based on the same period (30 vortex shedding cycles) used for the mesh and time-step refinement study in Table 1. This parameter shows the average force fluctuations at the two decks with respect to the average total value, and helps to understand the flow physics as a function of AoA. In Figure 15, the windward deck has quite stable RSD values within the investigated range of AoA, i.e., RSD values of 5% to 10%. For the cases of AoA = +6.3° and +8.1°, RSD values are less than 1%, see red line plot in Figure 10. The leeward deck has much larger RSD values within the investigated range of AoA as compared to the windward deck. The values show that the RSD values decrease as AoA increases. In Figure 16, positive lift forces on the windward deck can be observed when AoA ≥ 1.5°, whereas the lift force on the leeward deck is negative at all AoA. The smallest RSD magnitude of C_L on windward deck occurs at AoA = 6.3°, and the smallest RSD magnitude of C_L on leeward deck occurs at AoA = 8.1°, which indicate the change of the flow physics at these two AoA. The RSD of C_L on the windward deck at AoA = 0.1° and the RSD of C_L on the leeward deck at AoA = 6.3° show extremely large magnitude (i.e., 1627% and 977% respectively) due to the small value of $\overline{C_L}$ at these two AoA (see also Figure 8 and Figure 10).

4. Conclusions

The aerodynamic force coefficients and the Strouhal number for a preliminary design of the long-span bridge girder cross-section have been obtained numerically, using the k–$\omega$ SST turbulence model at Re$=$31,000. The objective of this study is to examine the validity of the 2D URANS simulations for the flow around a sharp-edged twin-box cross section for the quick and robust evaluation of the bridge design. The aerodynamic coefficients have been calculated and compared with the wind tunnel test results. The vortex formation mechanisms as a function of AoA have been studied. The relative contribution of each deck to the resultant force has also been investigated. The main conclusions are summarized as follows:

• The numerically predicted time-averaged force coefficients $\overline{C_D}$, $\overline{C_L}$ and $\overline{C_M}$ are in a good agreement with the wind tunnel experiment results. In particular, the drag coefficient has shown a good agreement with the experimental data at different AoA. The lift and moment coefficients show a good agreement with the experimental measurement in the low AoA range. A large reduction of $\overline{C_L}$ and $\overline{C_M}$ is observed at AoA = +10.2°. This indicates a premature stalling of the bridge decks simulated by the turbulence model, compared to the experimental observations in this high AoA region. The discrepancy may be due to the three-dimensional effects along the spanwise direction, which cannot be captured using the present 2D numerical model. Such high angles of attack are not expected during the normal bridge operation. Thus, the present 2D simulations are generally able to provide efficient and reliable assessment of the bridge girder aerodynamic performance under normal operating conditions.
• The flow structure shows a different pattern from AoA = +6.3°, as the vortices of the upstream deck become rather steady, while the vortices traveling along the upper surface of the downstream deck begin to merge towards a complete flow separation at even higher angles of attack. The flow pattern variations at AoA = +10.2°, −8.1° and −10.0° also have been discussed.
• Relative contributions of each deck to C_D and C_L varies with the AoA. This is also a good assessment tool, aiding the understanding of the flow physics and the screening purpose of the bridge design.

Overall, the present 2D URANS simulations with k–$\omega$ SST turbulence model give satisfactory results compared with the experimental results. The investigation of the flow pattern in the set of CFD studies can provide key information and insights to judge feasibility of aerodynamic performance of the bridge in a relatively cost-efficient manner.