Shape optimization of closed-box girder considering dynamic and aerodynamic effects on flutter: a CFD-enabled and Kriging surrogate-based strategy

Abstract

Shape optimization of single box girder allows the bridge to achieve the most favourable aerodynamic performance. Currently, aerodynamic optimization of bridge deck is usually achieved by traversing all potential schemes using wind tunnel tests or computational fluid dynamics (CFD), which are time-consuming, laborious, and restricted to few numbers of geometric parameter. This study proposed an aerodynamic optimization strategy using CFD technique, uniform design, Kriging surrogate model and hybrid infill sampling criteria to achieve the most favourable girder shape with the largest flutter wind speed. The effect of the variation of section shape on the dynamic characteristics of a real suspension bridge is considered. The Kriging surrogate model is then applied to approximately describe the aerodynamic parameters including static force coefficients and flutter derivatives as well as dynamic characteristics of the bridge. The cross-validation and samples infill are conducted to validate and improve the accuracy of the model. The optimal ranges, influence degree of each factor and their interactions on flutter performance are examined by statistical analysis. The contributions of aerodynamic parameters and dynamic characteristics to flutter performance are compared and discussed. The girder shape with the best flutter performance is obtained in updated surrogate model.

KEYWORDS:

• Streamlined closed-box girder
• CFD
• Shape Optimization
• Kriging surrogate model
• Uniform design
• Flutter

1. Introduction

The continuous advancements of construction technology, new materials, and design theories enable the increase of the span length of the cable-supported bridge to cross the wide straits and canyons. The wind-induced flutter critical velocity is getting lower, which is one of the great challenges for these bridges to survive in some extreme wind climates (Hu et al., 2023), such as strong typhoons (Fang et al., 2021a, 2021b). The design of long-span bridges must ensure that the flutter critical velocity is sufficiently higher than the extreme wind speed in a certain return period at the bridge site to avoid the occurrence of the devastating flutter instability (Montoya et al., 2018a; Fang et al., 2022). To further improve the flutter performance, the structural, mechanical or aerodynamic countermeasures are usually utilized. The structural measures are used to improve the stiffness of the bridge, such as spatial cables, central buckles, and wind-resistant cables, which are cost-consuming and inconvenient for design and operation. The mechanical measures, such as the tuned mass damper can also be applied to improve the flutter instability wind speed. But most mechanical measures are applicable to the motion in single degree of freedom (DoF) which are challenging to be used in two DoFs or three DoFs coupled flutter vibration. Moreover, flutter onset is insensitive to the damping of the structure, suggesting that the use of the dampers is also unsatisfactory (Abbas & Morgenthal, 2016). Comparatively, aerodynamic countermeasure is saving in cost and convenient in implementation by tailoring the deck configuration or installing some subsidiary components (Zhao et al., 2019). Although there are many passive aerodynamic countermeasures can be installed around the bridge deck to obviously enhance the flutter capacity of the bridge, such as the central stabilizer (Chen et al., 2006), the potential of the bridge girder to achieve an optimal aerodynamic performance by tailoring its geometric shape still deserve to be excavated. A tiny modification of the deck shape could significantly change the flutter critical wind speed since the flow separation and reattachment around the bluff body is very sensitive to its configuration (Wu et al., 2012; Siedziako & Øiseth, 2017).

The wind tunnel test is the widely used technique to investigate the flutter performance of the bridge deck with different shapes. Wang et al. (2009) studied the flutter performance of the section with different angles and widths of the wind fairing, showing that the wider and more sharpened wind fairing can strengthen the aerodynamic stability. Based on a large number of wind tunnel tests, Larsen and Wall (2012) recommended that the angle of inclined lower web should be less than 15° to suppress the vortex-induced vibration. Yang et al. (2018) found that the best flutter performance occurs at the angles of inclined lower web of 16° and 18° when the angles of the wind fairing are 50°and 60°, respectively. He et al. (2017) discovered that the effects of the aspect ratio, wind fairing angle and relative height of wind nose on the aerodynamic performance of the girder are negligible at small angles of wind attack, i.e. $|\alpha |$< 4°, but the flutter performance is very sensitive to these geometric parameters at $|\alpha |$> 4°. Although these wind tunnel tests provide many beneficial outcomes, they are confined to some specified geometrical parameters within a limited range. The interaction effects among different parameters on the flutter performance are not fully understood. The requirement of a large amount of time and cost as well as the manufacture of many sectional models make the wind tunnel tests inapplicable to solve a multi-parameter optimization problem. Alternatively, the continuous increase of computation capacity allows the computational fluid dynamics (CFD) technique to perform the aerodynamic shape optimization of the bridge deck during the preliminary design stage. Li et al. (2018) found that the wind fairing angle has insignificant effects on the flow field around the lower surface of streamlined box girder, but strikingly affects that around the upper surface. Haque et al. (2016a) studied the influence of fairing angle on aerodynamic response based on CFD simulations and found that the bottom plate slopes of 15°∼25° and top plate slopes of 40° or less lead to lower aerodynamic responses. Nevertheless, a too small bottom plate slope, say 10° would deteriorate the aerodynamical performance. Furthermore, Haque et al. (2016b) also explored the aerodynamic features of a pentagonal-shaped bridge deck by CFD using two-dimensional unsteady RANS (Reynolds-Averaged Navier-Stokes) simulation. The aerodynamic behaviour is clarified at high wind speeds by varying Reynolds number (2 × 10⁴–20 × 10⁴). It is found that the trailing-edge flow separation stops at the bottom plate slope of 11° with the Reynolds number of 3.9 × 10⁴, the flow separation may stop at the bottom plate slope greater than 11° with Re ≥10⁵.

Although CFD can make up for some shortcomings of the wind tunnel test, significant computing resources are still required, e.g. m^k full factor calculations for optimizing a box girder with m levels of k factors. For the preliminary design, what we need is to find an optimal section with the best flutter performance using as few computation resources as possible. The parameter domain with low flutter critical wind speeds should be identified first to reduce the computation costs in these regions. To achieve this, the surrogate model can be introduced to approximately describe the variation trend of the aerodynamic performance within the whole parameter domain using a small number of samples. In the field of aerospace and vehicle, the method of combining the surrogate model with CFD aerodynamic optimization design (Han et al., 2013) and structural optimization (Raponi et al., 2019) has been widely applied in high-speed trains (Xu et al., 2017), aircraft wing fuselage (Liu et al., 2022), and anti-icing system (Ridha et al., 2014), missile design (Vidanović et al., 2017) and unmanned combat aerial vehicles (Kaya et al., 2019).

The commonly used surrogate models include polynomial response surface method, Kriging method, support vector machine and so on. Among them, the Kriging method has a good approximation ability to nonlinear functions and allows a unique error estimation (Bernardini et al., 2015; Jeong et al., 2005), which will be used in this study. Besides, a reasonable sampling plan in the design domain is of great importance to the development of the surrogate model. The Latin hypercube sampling and uniform design are two widely used approaches. Montoya et al. (2018b) proved the feasibility and accuracy of the method combining CFD, Latin hypercube sampling and Kriging surrogate model to optimize the flutter performance of a closed-box girder with two geometrical parameters: width and depth of section. On this basis, Nieto et al. (2020) studied the flutter velocity response surface of a twin-box girder with three parameters of the gap distance, depth of section and the angle of wind fairing. Montoya et al. (2018a) also presented an approach to optimize the deck shape and cable sizes of a long-span cable-stayed bridge simultaneously considering the aeroelastic and structural constraints, including the cross-section area and prestressing force of each stay cable, the deck plates thickness and the width and depth of the box deck. Montoya et al. (2021) updated methodology to minimize the required amount of material of the bridge while meeting the performance and safety requirements in the structural and flutter design constraints. Sun et al. (2021) further introduced the ant colony genetic algorithm to optimize the flutter performance of a single box girder using Kriging surrogate model and uniform design. Considering the probabilistic flutter constraint, Kusano et al. (2020) performed a reliability-based optimization design for a long-span suspension bridge accounting for the effects of the shape and thickness of the girder. For the closed-box girder, above-mentioned studies mainly concern two or three geometrical shape parameters within a limited domain, which is unable to reveal the global characteristics of the flutter performance with respect to more parameters. And the degree of influence of each parameter is not well understood.

Moreover, the modification of the deck shape inevitably changes its cross-section property as well as the dynamic characteristics of the bridge, and finally affect the flutter critical velocity. Zhan and Liao (2019) found that the depth of deck has insignificant effects on the vertical bending frequency but obviously increases the torsional frequency. Nieto et al. (2011) found that the natural frequencies and mode shapes are mainly influenced by the cross-section area or the deck mass. Moreover, sensitivities of the flutter velocity with respect to cross-section area and torsional inertia are positive, to the bending inertias are negative. Wang et al. (2014) also found that the increase of torsional stiffness would increase the torsional frequency and improve the flutter stability by conducting the sensitivity analysis. But the vertical stiffness has neglectable effects on flutter performance. Abbas and Morgenthal (2016) implemented sensitivity analysis on the uncertainty of flutter derivatives and structural parameters, showing that the frequency ratio has a comparable effect on flutter to that of key flutter derivatives. These studies have revealed that the sensitivity of the flutter critical velocity to the dynamic characteristics could be equally important to the aerodynamic configuration.

The main objective of this study is to establish a shape optimization process of the single box girder considering aerodynamic and dynamic effects on the flutter performance of a long-span suspension bridge. In Section 2, the theoretical background of bridge flutter, uniform design and surrogate models is introduced. Section 3 describes the optimization process and selection scheme of design factors, design domain and samples. In Section 4, the CFD numerical simulation and independence test are performed before introducing the finite element model of a real bridge. The surrogate models of static force coefficients and their slopes, flutter derivatives and fundamental frequencies are then presented in Section 5, respectively. The accuracy of surrogate model is improved and tested by cross-validation and samples update. The optimal ranges of five dimensionless parameters for the closed-box girder are recommend. The statistical analysis of the influence trend and importance of each factor are also implemented. The contributions of dynamic characteristics and aerodynamic geometrical shape to flutter performance are compared and discussed. Finally, the optimal and suboptimal cross-sectional shape parameters with maximum flutter critical velocities are determined.

2 Theoretical background

2.1 Linear flutter model of bridge section

To describe the flutter motion of the bridge deck, the linear self-excited aerodynamic force modelled as the function of the motion displacement and velocity of the bridge deck is commonly used. The flutter derivatives are introduced to connect the motion and surface force of the bridge deck. For a streamlined box girder, the contribution of the drag component to the self-excited force is insignificant and negligible, and the corresponding lateral degree of freedom (DoF) is ignored. The equations of motion related to vertical and torsional DoFs are given by (1) $\begin{array}{rl}& \ddot{h}+2{\xi }_{h0}{\omega }_{h0}\dot{h}+{\omega }_{h0}^{2}h\\ & =\frac{\rho B{U}^2}{m_h}\left(K{H}_1^{\ast }\frac{\dot{h}}{U}+K{H}_2^{\ast }\frac{B\dot{\alpha }}{U}+K^2{H}_3^{\ast }\alpha +K^2{H}_4^{\ast }\frac{h}{B}\right)\end{array}$(1) (2) $\begin{array}{rl}& \ddot{\alpha }+2{\xi }_{\alpha 0}{\omega }_{\alpha 0}\dot{\alpha }+{\omega }_{\alpha 0}^2\alpha \\ & =\frac{\rho B^2{U}^2}{I_m}\left(K{A}_1^{\ast }\frac{\dot{h}}{U}+K{A}_2^{\ast }\frac{B\dot{\alpha }}{U}+K^2{A}_3^{\ast }\alpha +K^2{A}_4^{\ast }\frac{h}{B}\right)\end{array}$(2) where U is wind velocity, B is the width of section, $K=\omega B/f$ is reduced frequency, f is vibration frequency. $h$ and $\dot{h}$ are vertical bending displacement and velocity, $\alpha$ and $\dot{\alpha }$ are torsional displacement and velocity, $H_i^{\ast }\text{ and }{A}_i^{\ast }$ (i = 1, 2, 3, 4) are the flutter derivatives. ${\omega }_{h0}=2\pi f_{h0}$ and ${\omega }_{\alpha 0}=2\pi f_{\alpha 0}$ are the circular frequencies of vertical and torsional modes, respectively, ${\xi }_{h0}$ and ${\xi }_{\alpha 0}$ are the corresponding modal damping ratios. In the flutter critical state, the vibration of vertical and torsional DoFs can be formulated as (3) $h=h_0{e}^{i\omega t}\alpha ={\alpha }_0{e}^{i\left(\omega t+\theta \right)}$(3) where $h_0$ and ${\alpha }_0$ are initial amplitudes, $\theta$ is the vibration phase difference between vertical and torsional DoFs. Substituting Eq. (3) into Eqs. (1) and (2), the vibration frequency and damping ratio of two vibration modes at different velocities can be iteratively obtained. For a streamlined bridge section, the frequency and damping ratio related to torsional mode will dominate the occurrence of the flutter, which can be formulated as: (4) $\begin{array}{rl}{\omega }_{\alpha }& ={\omega }_{\alpha 0}/\{1+\frac{\rho B^4{A}_3^{\ast }}{I_m}+\frac{{\rho }^2{B}^6}{I_m\cdot m_h}\\ & \times {\mathrm{\Omega }}_{h\alpha }\left[\right(A_4^{\ast }{H}_2^{\ast }+A_1^{\ast }{H}_3^{\ast })sin{\theta }_{h\alpha }\\ & +{\left(A_4^{\ast }{H}_3^{\ast }-A_1^{\ast }{H}_2^{\ast }\right)cos{\theta }_{h\alpha }\left]\phantom{\frac{{\rho }^2{B}^6}{I_m\cdot m_h}}\right\}}^{1/2}\end{array}$(4) (5) $\begin{array}{rl}{\xi }_{\alpha }& =\frac{{\xi }_{\alpha 0}{\omega }_{\alpha 0}}{{\omega }_{\alpha }}-\frac{\rho B^4{A}_2^{\ast }}{2I_m}-\frac{{\rho }^2{B}^6}{2I_m\cdot m_h}\\ & \times {\mathrm{\Omega }}_{h\alpha }\left[\right(A_1^{\ast }{H}_2^{\ast }-A_4^{\ast }{H}_3^{\ast })sin{\theta }_{h\alpha }\\ & +\left(A_4^{\ast }{H}_2^{\ast }+A_1^{\ast }{H}_3^{\ast }\right)cos{\theta }_{h\alpha }]\end{array}$(5) where ${\mathrm{\Omega }}_{h\alpha }$ and ${\theta }_{h\alpha }$ can be obtained from (6) $\begin{array}{rl}{\mathrm{\Omega }}_{h\alpha }& =\frac{{\omega }_{\alpha }^2}{\sqrt{{\left({\omega }_h^{{}^{\mathrm{\prime }}2}-{\omega }_{\alpha }^2\right)}^2+{\left(2{\xi }_h^{\mathrm{\prime }}{\omega }_{\alpha }^2\right)}^2}}\end{array}$(6) (7) $\begin{array}{rl}{\theta }_{h\alpha }& =\arctan \frac{2{\xi }_h^{\mathrm{\prime }}{\omega }_{\alpha }^2}{{\omega }_h^{{}^{\mathrm{\prime }}2}-{\omega }_{\alpha }^2}+\pi \end{array}$(7) where ${\theta }_{h\alpha }$ is phase difference of vertical and torsional vibration displacements. For long-span bridges, ${\omega }_{\alpha }$ is usually larger than ${\omega }_h$, and $\frac{\pi }{2}<{\theta }_{h\alpha }<\pi$. With the increase of wind velocity, ${\omega }_{\alpha }$ and ${\xi }_{\alpha }$ are iteratively solved by Eqs. (4) and (5). The flutter critical velocity can be solved when ${\xi }_{\alpha }=0$.

2.2 Flutter derivatives

The estimation of flutter derivatives for the bridge deck provides basic inputs to predict the flutter onset. Wind tunnel test or CFD simulations using free decay vibration or forced vibration methods is widely employed. For the convenience of operation and reduction of computation amount, the quasi-steady assumption is adopted in this study. The quasi-steady theory (Scanlan, 1988) is a modelling scheme for wind-induced loads on bridge decks which takes into account deck displacement and wind turbulence properties. It relates the displacements of the bridge to the relative wind velocity component and consists of a two-dimensional framework based on the ‘strip theory’ (Davenport, 1962). Therefore, the model is adequate when the time required for the flow to travel far enough downstream across the bridge deck is much lower than the time required for the structure to respond to changes in the surrounding flow. So, this assumption is often used for high reduced wind velocity and streamlined deck sections (Tubino, 2005). The flutter derivatives can be approximately estimated using quasi-steady theory: (8) $\begin{array}{rl}{H}_1^{\ast }& =-\frac{C_{L,0^{\circ }}^{\mathrm{\prime }}}{2K}{H}_2^{\ast }\approx -\frac{C_{L,0^{\circ }}^{\mathrm{\prime }}}{2K}{\mu }_H{H}_3^{\ast }\approx -\frac{C_{L,0^{\circ }}^{\mathrm{\prime }}}{2K^2}{H}_4^{\ast }=0\end{array}$(8) (9) $\begin{array}{rl}{A}_1^{\ast }& =\frac{C_{M,0^{\circ }}^{\mathrm{\prime }}}{2K}{A}_2^{\ast }=\frac{C_{M,0^{\circ }}^{\mathrm{\prime }}}{K}{\mu }_A{A}_3^{\ast }\approx \frac{C_{M,0^{\circ }}^{\mathrm{\prime }}}{2K^2}{A}_4^{\ast }=0\end{array}$(9) (10) $\begin{array}{rl}K& =\frac{B\omega }{U}{\mu }_A\approx \frac{A_3^{\ast }}{H_3^{\ast }}{\mu }_H\approx \frac{A_1^{\ast }}{H_1^{\ast }}\end{array}$(10) where $C_{D,0^{\circ }}$, $C_{L,0^{\circ }}$, $C_{M,0^{\circ }}$, $C_{L,0^{\circ }}^{\mathrm{\prime }}$ and $C_{M,0^{\circ }}^{\mathrm{\prime }}$ are the aerodynamic coefficients and their first-order derivatives with respect to the angle of attack at α = 0°. It is noteworthy that the drag, lift and moment coefficient in these equations are normalized based on the width of section.

2.3 Uniform design

The uniform design was proposed by Wang and Fang (1981), which considers the index of ‘uniform dispersion’ and ‘space filling’. The design domain space information can be obtained to the maximum extent by selecting the fewest samples. Uniform design tables (UDT) are divided into two types, including ${\text{UD}}_{\text{n}}^{\ast }$(q^m) and ${\text{UD}}_{\text{n}}^{\ast }$ (q^m), where * means better uniformity and should be preferable, UD refers to the uniform design, n represents the number of tests, q refers to the level of each factor, and m refers to the number of factors. Each UDT should be attached a usage table that explains how to select columns for the UDT.

The UDT can be formed by Latin method, the threshold-accepting method and integer programming problem method (Fang et al., 2018). The good lattice point method (GLP) (Dai & Wang, 2009) is considered to be the most widespread method which is introduced here. Firstly, find an integer h less than the sample number n with the greatest common divisor between them of 1. All positive integers meet the previous conditions will give us a generator vector h = (h₁, h₂, … , h_m), where m is the maximum column number of a UDT can have. Then the j column of UDT is generated according to the following equation: (11) $u_{ij}=i{h}_j\left[\text{mod}n\right]$(11) where [mod n] represents the congruence operator, u_ij refers to the corresponding element of row i and column j of UDT. If ih_j exceeds n, then an appropriate multiple of n is subtracted by ih_j so that the difference falls in [1, n]. And u_ij can be generated recursively: (12) $\begin{array}{rl}{u}_{1j}& =h_j\end{array}$(12) (13) $\begin{array}{rl}{u}_{i+1,j}& =\{\begin{array}{l}{u}_{ij}+h_j\left(u_{ij}+h_j\right)\le n\\ u_{ij}+h_j-n\left(u_{ij}+h_j\right)>n\end{array}i=1,2,\cdots ,n-1\end{array}$(13) The uniformity of UDT is measured by the discrepancy G. The smaller the value of discrepancy G, the better the uniformity of the sample. Transform n trial points into n points in [0, 1] = C^m: x1, … , xn. The uniformity of the original n test points is equivalent to measuring the uniformity of x1, … , xn in Cm. G is formulated as (14) $G\left({\mathit{x}}_1,\dots ,{\mathit{x}}_n\right)=\text{sup}\left\{\mathit{x}\in C^m\right\}\left|\frac{n_{\mathit{x}}}{n}-v\left(\mathit{x}\right)\right|$(14) where ${\text{x}}_{\text{1}},\dots ,{\text{x}}_{\text{n}}$ stand for n points in ${\text{C}}^{\text{m}}$, for any vector x = $\text{(}{\text{x}}_{\text{1}},\dots ,{\text{x}}_{\text{n}})\in {\text{C}}^{\text{m}},$, v(x) = ${\text{x}}_{\text{1}},\dots ,{\text{x}}_{\text{n}}$ is the volume of rectangle $\left\{0,\mathit{x}\right\},n_{\mathit{x}}$, is the number of points falling into $\left\{0,\mathit{x}\right\}$ in ${\text{x}}_{\text{1}},\dots ,{\text{x}}_{\text{n}}$, and sup means the supremum.

2.4. Surrogate model

The basic idea of the surrogate model is to replace the complex high-precision analysis model with a simple approximation function at certain accuracy conditions. The response value of the known point (sample) is used to predict that of the unknown point. Among the existing surrogate models, Kriging model developed from geostatistics is a commonly used method, which was proposed by Krige (1951). Kriging model is an unbiased estimation model (Lophaven et al., 2002), which can not only provide the estimated value of the unknown function, but also provide the error estimation of the estimated value. This is the indigenous feature of Kriging model different from other surrogate models. It has been applied to solve some complex problems (Rahman et al., 2010; Roy et al., 2008), which can achieve great optimization efficiency and ensure the optimization quality.

In Kriging model, the prediction response $\hat{y}\left(x\right)$ and the Root Mean Square Error (RMSE) $\phi \text{(x)}$ of the prediction can be expressed as (15) $\begin{array}{rl}\hat{y}\left(\mathit{x}\right)& =f(\mathit{x}{)}^T\text{β}+z\left(\mathit{x}\right)\end{array}$(15) (16) $\begin{array}{rl}{\phi }^2\left(\mathit{x}\right)& ={\mathit{\sigma }}^2\left(1+{\mathit{c}}^T\mathit{R}\mathit{c}-2{\mathit{c}}^T\mathit{r}\right)\end{array}$(16) where $\text{β}$ is a linear regression function and $f\text{(x)}$ is a polynomial function of variable x, which provide a global approximation for the establishment of the model. ${\sigma }^{\text{2}}$ is a proportional factor, which is called process variance. R is the correlation matrix of variable x. r represents the correlation between the sample points and the prediction points. c offers the weight of the linear combination in the Kriging surrogate model and should minimize the RMSE. z(x) is a random process with non-zero covariance, which follows the normal distribution N(0, ${\sigma }^2$) and offers a local approximations of model deviations. The covariance matrix of z(x) can be written as (17) $\begin{array}{rl}E\left[z\left(\mathit{w}\right)z\left(\mathit{x}\right)\right]& ={\mathit{\sigma }}^2\mathcal{R}\left(\mathit{\theta },\mathit{w},\mathit{x}\right)\end{array}$(17) (18) $\begin{array}{rl}\mathcal{R}\left(\mathit{\theta },\mathit{w},\mathit{x}\right)& =\prod _{j=1}^n{\mathcal{R}}_j\left(\mathit{\theta },w_j-x_j\right)\end{array}$(18) (19) $\begin{array}{rl}{\mathcal{R}}_j\left(\mathit{\theta },w_j-x_j\right)& =\text{exp}\left(-{\theta }_j{|w_j-x_j|}^2\right)\end{array}$(19) where $\mathcal{R}(\theta ,\text{w},\mathrm{x})$ is the spatial correlation function of any two-sample points w and x.

The parameter $\mathit{\theta }$ in the correlation function is unknown to be estimated, which represents the correlation of samples in various spatial dimensions. n is the dimension of the input parameter. Gaussian kernel function commonly used in ${\mathcal{R}}_{\text{j}}(\theta ,{\text{w}}_{\text{j}}-{\text{x}}_{\text{j}})$ is employed in this study, which has the best calculation effect (Giunta et al., 1997).

2.5. Hybrid infill sampling criteria

The accuracy of the first generation of the surrogate model is usually inadequate due to insufficient samples or nonlinearity of the problem. To further increase the accuracy, surrogate model is usually updated by introducing more samples. The hybrid infill sampling criteria, including the expected improvement (EI), probability of improvement (PI), minimize prediction (MP) and root mean square error (RMSE) are employed. The EI criterion accounts for the effects of the predicted value of the function and the error of the prediction achieved by the Kriging model. The EI value of the Kriging model at any unknown point x can be formulated as: (20) $\begin{array}{rl}EI\left(\mathit{x}\right)& =\left[\hat{y}\left(\mathit{x}\right)-y_{max}\right]\mathrm{\Phi }\left(\mathit{Z}\right)+\phi \left(\mathit{x}\right)\varphi \left(\mathit{Z}\right)\end{array}$(20) (21) $\begin{array}{rl}\mathit{Z}& =\{\begin{array}{ll}{\displaystyle \frac{\hat{y}\left(\mathit{x}\right)-y_{max}}{\phi \left(\mathit{x}\right)}}& if\phi \left(\mathit{x}\right)>0\\ 0& if\phi \left(\mathit{x}\right)=0\end{array}\end{array}$(21) where $\mathrm{\Phi }$(·) and $\varphi$(·) indicate the CDF and PDF of the standard normal distribution, respectively. ymax is the maximum value in the known samples.

The PI criterion assumes that you want a larger value (T) than ${\hat{y}}_{max}$, then the probability of greater than T is as follows: (22) $\begin{array}{rl}T& ={\hat{y}}_{max}+\kappa |{\hat{y}}_{max}|\kappa \ge 0\end{array}$(22) (23) $\begin{array}{rl}PI\left(\mathit{x}\right)& =\mathrm{\Phi }\left(\frac{\hat{y}\left(\mathit{x}\right)-T}{\phi \left(\mathit{x}\right)}\right)\end{array}$(23) where $\kappa$ = 0.1, i.e. a solution expected to be 10% larger than the current true optimal value.

As compared with the widely used criterion of minimizing prediction, the criterion of maximizing the predicted objective function (also called MP) is used: (24) $MP=\hat{y}(\mathit{x}{)}_{max}$(24) The Root Mean Square Error (RMSE) criterion adds samples where the model error is the largest. (25) $RMSE\left(\mathit{x}\right)=\phi \left(\mathit{x}\right)$(25)

The above four criteria are summarized in Table 1. EI and PI are global optimization infill sampling criteria, which are able to find the global optimal value if sufficient samples are infilled (Törn & Žilinskas, 1989; Locatelli, 1997). However, EI and PI near the true global optimum could be small after infilling many samples, resulting in that the subsequent infill samples are mostly in the sparse region or other local optimal regions, thereby reducing the convergence speed and accuracy. MP is a local optimization infill sampling criterion, but it probably falls into the local optimal value and fails to find the global optimum. To balance the convergence speed and accuracy of the prediction, the hybrid infill sampling criteria is used in this paper. EI and PI criteria are used to find global optimal value in design domain. MP is used to more accurately locate the extreme value point near the current optimal value. RMSE is used to improve the global accuracy of surrogate models and provides more accurate $\hat{y}\left(\mathit{x}\right)$ and $\phi \left(\mathit{x}\right)$ for other infill sampling criteria. The positions of max(EI(x)), max(PI(x)), MP and max(RMSE(x)) are infilled with new samples.

Table 1. Characteristics of 4 types of infill sampling criteria.

Criterion	Application	Virtue	Defect
EI	Global optimization	Find the global optimum if infilling sufficient samples	Low convergence speed and poor accuracy
PI	Global optimization	Similar to EI (different optimization algorithm)	Similar to EI
MP	Local optimization	Fast convergence speed and more accuracy	Falls into the local optimal value
RMSE	Improve the global accuracy	Improve the overall accuracy	No significant effect on the optimal value

3 Optimization preparation

3.1. Optimization process

The aerodynamic optimization of closed-box girder refers to the parameterization of geometrical shape such as width, and depth of the section, the wind fairing, inclined upper and lower webs, which are called design parameters (factors). The appropriate variation range for each parameter is selected according to the parameter constraints in the actual design, which is called the design domain. The interval size of each parameter in the design domain needs to be able to capture the target value, and the interval is determined, namely, the level number is determined. The optimization goal in this study is to obtain a closed-box girder with the maximum flutter critical wind speed.

The optimization process is summarized in Figure 1: (a) Preliminary preparation: Selecting design parameters, design domain and levels before determining the experimental scheme or experimental points using uniform design. (b) Numerical calculation: (1) Developing the aerodynamic parameters, including CD, CL, CM, ${\text{C}}_{\text{D}}^{\mathrm{\prime }}{\text{, C}}_{\text{L}}^{\mathrm{\prime }}{\text{, C}}_{\text{M}}^{\mathrm{\prime }},{\text{H}}_{\text{i}}^{\ast }{\text{ and A}}_{\text{i}}^{\ast }\text{ (i = 1,2,3,4)}$ using the results from the CFD simulation. (2) Developing the dynamic characteristics, including the vertical bending frequency fh and torsional frequency fα of first-order symmetric and anti-symmetric modes as well as their equivalent mass m and equivalent inertia moment Im per unit length using the results from the finite element analysis. (c) Construction of Kriging surrogate model for flutter critical velocity: Calculating the flutter critical velocities based on the parameters predicted by the Kriging approximation method. (d) Cross validation: More sample points will be updated in the parameters domain using hybrid infill sampling criteria before reaching a certain accuracy. (e) Statistical analysis: The main effect of each factor, the influence of each factor and their interaction on flutter performance are analyzed. The contributions of aerodynamic parameters and dynamic characteristics to flutter performance are discussed. (f) Optimal Section: The configuration parameters with maximum critical wind speed are then determined to obtain the optimal section.

Figure 1. Aerodynamic optimization process of closed-box girder.

It is noteworthy that the bare cross-section without any auxiliary facilities is employed in this study, which is mainly utilized in the preliminary design stages. With the continuous increase of computation capacity, we believe that the aerodynamic shape optimization of the real bridge deck equipped with all auxiliary facilities can be achieved in the future. But the methodology proposed in this study is also applicable.

3.2 Design parameters and domain

Many studies (Wang et al., 2009; Larsen & Wall, 2012; He et al., 2017; Yang et al., 2018; Zhao et al., 2019) have revealed that the angles of wind fairing and inclination of the lower web have significant effects on the aerodynamic performance of the streamlined closed-box section. However, these two angle parameters are not independent, which are inappropriately used as design factors. Alternatively, the coordinates of the wind faring corner point and the intersection point between the inclined lower web and bottom panel are utilized to describe the configuration the bridge deck. The width of the top deck panel remains unchanged to meet the requirement of the traffic volume.

A closed-box girder section with eight-lane dual carriageway is employed as the benchmark. The width of top plate $B^{\mathrm{\prime }}$ = 40.5 m and transverse slope of 2.4%. The vertical and horizontal coordinates of the wind fairing and inclined lower web points (points A and B in Figure 2(b)) are selected as the four optimization parameters or design factors. The horizontal coordinates of points A and B can be described using the width of section and bottom plate. Similarly, the vertical coordinates of these two points are expressed as the vertical distances away from the central point of the bridge deck, which are defined as H and D, respectively.

Figure 2. Variation range of aerodynamic configuration.

Specifically, Figure 2 visualizes the size of the design domain. The design domain of point A is within a yellow point line box with the range of horizontal coordinate B1 from 41.8 m to 58.6 m and the vertical coordinate H from 1.05 m to 3.15 m. Point B refers to the intersection point between the inclined lower web and bottom deck, whose design domain is a blue dotted box. The range of horizontal coordinate for point B is from 24.5 m to 35.7 m, while the vertical coordinate D varies from 2.88 m to 6.24 m. Within the design domain, the width-depth ratio B/D = 6.70−20.35. The top θ1 and bottom θ2 represent to the angles between the inclined sides of the deck to the horizontal axis, and θT is the complete angle of the wind fairing, i.e. θT = θ1+θ2. The values of three angles in the design domain are shown in Table 3, where θ1 = 5.03°∼71.57°, θ2 = −5.06°∼59.56° and θT = 12.35°∼122.74°.

3.3. Experimental design

After the design domain is determined, samples are selected using the uniform design method described in Section 2.3. Since there are four design parameters in this study, the levels of each factor should be as dense as possible to ensure that the discrepancy G is sufficiently small and the resolution of the flutter velocity is high enough to capture the maximum velocity region. The ${\text{UD}}_{\text{29}}^{\ast }$(296) UDT is utilized, and its usage table is shown in Table 2, where column No.1, 3, 4, and 5 are selected with discrepancy G = 0.105. Especially, the s refers number of design factors, namely, number of columns selected of the UDT. The initial experimental design scheme is shown in Figure 3, in which 68 groups of shape schemes (black circles) are selected for uniform design (Note: 29 groups were obtained by looking up the UDT, and the remaining 39 groups expanded the range of three parameters (H, D, B₂) by half based on uniformity). The values of parameters for each sampling plan are shown in Table 3. The initial experimental design scheme is shown in Figure 3, in which 68 groups of shape schemes (black circles) are selected for uniform design. The values of parameters for each sampling plan are shown in Table 3.

Figure 3. Uniform experimental design: sample distribution.

Table 2. Usage table of ${\text{UD}}_{\text{29}}^{\ast }$(29⁶) UDT.

s	Column No.					discrepancy G
2	1	4				0.052
3	1	3	4			0.091
4	1	3	4	5		0.105
5	2	3	4	5	6	0.173

Table 3. Basic parameters of numerical model.

No.	B2	B1	H	D	B1/D	${\theta }_1{/}^{\circ }$	${\theta }_2{/}^{\circ }$	${\theta }_T{/}^{\circ }$	No.	B2	B1	H	D	B1/D	${\theta }_1{/}^{\circ }$	${\theta }_2{/}^{\circ }$	${\theta }_T{/}^{\circ }$
1	24.5	48.2	1.95	4.64	10.39	20.64	12.79	33.43	35	28.9	51	1.9	5.68	8.98	14.93	18.88	33.82
2	24.9	43	1.4	4.08	10.54	35.75	16.50	52.25	36	30.5	42.2	1.2	5.84	7.23	39.47	38.42	77.89
3	25.3	49.8	2.35	3.52	14.15	21.70	5.46	27.15	37	30.9	49	2.15	5.28	9.28	21.22	19.08	40.30
4	25.7	44.6	1.8	2.96	15.07	32.38	7.00	39.38	38	32.1	45.4	2	6	7.57	31.48	31.03	62.50
5	26.1	51.4	1.25	4.8	10.71	7.84	15.68	23.51	39	32.5	52.2	1.45	5.44	9.60	9.22	22.05	31.28
6	26.5	46.2	2.2	4.24	10.90	30.82	11.70	42.52	40	33.7	48.6	1.3	6.16	7.89	11.17	33.12	44.29
7	26.9	53	1.65	3.68	14.40	10.43	8.84	19.27	41	34.1	43.4	2.25	5.6	7.75	50.36	35.77	86.13
8	27.3	47.8	1.1	3.12	15.32	9.33	11.15	20.48	42	35.7	46.6	1.55	5.76	8.09	19.00	37.69	56.68
9	27.7	42.6	2.05	4.96	8.59	55.89	21.34	77.22	43	27.7	54.6	2.05	4.96	11.01	12.40	12.21	24.61
10	28.1	49.4	1.5	4.4	11.23	12.67	15.23	27.90	44	28.5	56.2	2.45	6.24	9.01	13.95	15.30	29.25
11	28.5	44.2	2.45	3.84	11.51	46.51	10.04	56.55	45	29.3	57.8	1.35	5.12	11.29	5.61	14.82	20.43
12	28.9	51	1.9	3.28	15.55	14.93	7.12	22.05	46	30.5	54.2	1.2	5.84	9.28	5.83	21.38	27.22
13	29.3	45.8	1.35	5.12	8.95	17.78	24.56	42.34	47	31.3	55.8	1.6	4.72	11.82	8.18	14.29	22.47
14	29.7	52.6	2.3	4.56	11.54	16.57	11.17	27.73	48	32.1	57.4	2	3.6	15.94	10.07	7.21	17.27
15	30.1	47.4	1.75	4	11.85	19.92	14.58	34.50	49	33.3	53.8	1.85	4.32	12.45	11.48	13.55	25.02
16	30.5	42.2	1.2	3.44	12.27	39.47	20.95	60.42	50	34.1	55.4	2.25	5.6	9.89	13.22	17.46	30.68
17	30.9	49	2.15	2.88	17.01	21.22	4.61	25.83	51	34.9	57	1.15	4.48	12.72	4.50	16.77	21.28
18	31.3	43.8	1.6	4.72	9.28	33.69	26.53	60.22	52	35.7	58.6	1.55	3.36	17.44	6.62	8.98	15.60
19	31.7	50.6	1.05	4.16	12.16	6.22	18.22	24.43	53	24.9	55	1.4	4.08	13.48	7.08	10.10	17.17
20	32.1	45.4	2	3.6	12.61	31.48	13.53	45.01	54	25.7	56.6	1.8	5.36	10.56	9.17	12.98	22.15
21	32.5	52.2	1.45	3.04	17.17	9.22	9.17	18.39	55	26.5	58.2	2.2	4.24	13.73	10.87	7.33	18.21
22	32.9	47	2.4	4.88	9.63	30.31	19.38	49.69	56	24.9	55	2.9	4.08	13.48	18.32	4.48	22.80
23	33.3	41.8	1.85	4.32	9.68	64.29	30.16	94.45	57	26.1	51.4	2.75	4.8	10.71	22.43	9.21	31.64
24	33.7	48.6	1.3	3.76	12.93	11.17	18.27	29.45	58	26.9	53	3.15	3.68	14.40	22.98	2.33	25.30
25	34.1	43.4	2.25	3.2	13.56	50.36	11.55	61.90	59	27.3	47.8	2.6	5.52	8.66	29.91	15.90	45.81
26	34.5	50.2	1.7	5.04	9.96	13.90	23.05	36.95	60	28.1	49.4	3	4.4	11.23	29.33	7.49	36.82
27	34.9	45	1.15	4.48	10.04	16.11	33.40	49.51	61	29.3	45.8	2.85	5.12	8.95	41.57	15.38	56.95
28	35.3	51.8	2.1	3.92	13.21	15.81	12.44	28.25	62	30.5	42.2	2.7	5.84	7.23	68.88	28.22	97.10
29	35.7	46.6	1.55	3.36	13.87	19.00	18.37	37.37	63	31.3	55.8	3.1	4.72	11.82	18.77	7.53	26.30
30	25.3	49.8	2.35	5.92	8.41	21.70	16.25	37.94	64	31.7	50.6	2.55	4.16	12.16	22.09	9.67	31.76
31	25.7	44.6	1.8	5.36	8.32	32.38	20.64	53.02	65	32.5	52.2	2.95	5.44	9.60	22.72	14.19	36.91
32	26.9	53	1.65	6.08	8.72	10.43	18.75	29.18	66	33.7	48.6	2.8	6.16	7.89	29.59	24.28	53.87
33	27.3	47.8	1.1	5.52	8.66	9.33	23.33	32.66	67	34.9	45	2.65	4.48	10.04	43.70	19.92	63.62
34	28.5	44.2	2.45	6.24	7.08	46.51	25.77	72.28	68	35.7	46.6	3.05	5.76	8.09	39.90	26.44	66.34

4 Numerical simulation

4.1 CFD numerical scheme

The POINTWISE 18.0 R1 and FLUENT 2021 R1 are used to mesh the grid and perform the calculation, respectively during the CFD numerical simulation. The sectional rigid model adopts a scale ratio of 1:80 and Reynolds number Re = 10⁵. The calculation domain is shown in Figure 4 (B and D are the width and depth of the section). The numerical calculation domain is 15B in depth and 37.5B in width. The centre of the streamlined closed-box girder section is 15B away from the inlet, 22.5B from the exit, and 7.5B from the upper and lower boundaries. The left inlet boundary of the computational domain is set as the velocity inlet 10 m/s. The right outlet boundary is set as the pressure outlet. The upper and lower boundaries are set as symmetric boundaries.

Figure 4. Design scheme of computing domain.

In order to ensure the reasonable density of the mesh as much as possible, the mesh block technology is used to divide the entire calculation domain into five zones. Structured meshes are applied to zone 1, zone 2, and zone 4 and unstructured meshes are used in zone3 and zone5. Forty-layer of structured meshes are used for the near wall surface with the bottom layer height of 5 × 10⁻⁵ m and the growth rate of 1.02–1.03. The maximum Y⁺ of all models is less than 2. The mesh number are from 240,000–250,000. The number of nodes in the bridge section is at range of 1294-1454. Figure 5 shows some details of mesh meshing by zooming in the mesh around the bridge section.

Figure 5. Partial meshing diagram.

The pressure-based SIMPLE solver, the spatial discretization of gradient and second-order interpolation and the second-order implicit time discretization are employed. The turbulence model adopts 2D k-ω SST. Application examples of this turbulence model can be found in Villalpando et al. (2012). The time step is 6 × 10⁻⁵s to ensure that the maximum Courant number of all models is less than 2.

4.2 CFD independence test and wind tunnel test validation

It is a necessary step to verify that the numerical simulation results are unaffected by mesh and time discretization. The tests of independence for mesh size and time step are performed, respectively. During the test, one of two parameters, i.e. the size of mesh and the time increment, is fixed by varying the other parameter to obtain an appropriate value. Three numbers of grid, i.e. roughly about 150,000, 250,000 and 340,000 are selected and the time increment Δt = 6 × 10⁻⁵s. The results are shown in Table 4. It can be seen that the 150,000-grid model has lower estimations of $C_D^{\mathrm{\prime }}$ and $C_L^{\mathrm{\prime }}$ compared with other two grid models. The rest of parameters are almost identical among these three grid numbers. In order to balance the computation efficiency and accuracy, the mesh number of 240,000 around is utilized.

Table 4. Force coefficients of different numbers of grid.

Total cells	Deck cells	Mesh layers	Maximum Courant number	$y_{max}^{+}$	$C_{D,0^{\circ }}$	$C_{L,0^{\circ }}$	$C_{M,0^{\circ }}$	$C_D^{\mathrm{\prime }}$	$C_L^{\mathrm{\prime }}$	$C_M^{\mathrm{\prime }}$
154877	1654	30	0.97	1.55	0.028	−0.066	0.009	0.020	5.181	1.458
249527	1654	40	0.98	1.55	0.028	−0.066	0.009	0.018	5.105	1.456
342585	1994	40	1.18	1.57	0.028	−0.065	0.009	0.018	5.114	1.457

Figure 6. TJ-7 wind tunnel.

Figure 7. Cross section of bridge deck.

Figure 8. Wind tunnel test of sectional model.

Similarly, three time increments are selected for the test of time increment independence, namely, Δt = 6 × 10⁻⁵s, Δt = 1.2 × 10⁻⁴s, and Δt = 3 × 10⁻⁴s. As shown in Table 5, slight differences of $C_D^{\mathrm{\prime }}$ and $C_L^{\mathrm{\prime }}$ can be found between Δt = 6 × 10⁻⁵s and other two time steps. In addition, the maximum Courant number for Δt = 6 × 10⁻⁵s is less than 1. Therefore, the time step of Δt = 6 × 10⁻⁵s is employed.

Table 5. Force coefficients of different time step.

Time step	Maximum Courant Number	$y_{max}^{+}$	$C_{D,0^{\circ }}$	$C_{L,0^{\circ }}$	$C_{M,0^{\circ }}$	$C_D^{\mathrm{\prime }}$	$C_L^{\mathrm{\prime }}$	$C_M^{\mathrm{\prime }}$
$6\times {10}^{-5}s$	0.98	1.55	0.028	−0.067	0.009	0.018	5.105	1.456
$1.2\times {10}^{-4}s$	1.95	1.55	0.028	−0.067	0.009	0.017	5.106	1.456
$3\times {10}^{-4}s$	4.73	1.55	0.028	−0.067	0.009	0.017	5.107	1.456

It is noteworthy that the use of POINTWISE is able to record the operation steps of meshing the grids, which allows the quick meshing of different bridge deck sections using their coordinates information. The batch meshing for all decks of interest can be achieved before submitting to the high-performance computers to conduct the parallel computation.

To verify the accuracy of the CFD simulation, wind tunnel tests are also conducted in the TJ-7 wind tunnel. TJ-7 is a straight flow air-blowing wind tunnel with designed wind speed of 0–15 m/s. The cross-sectional size is 0.65 m × 1.2 m. The maximum downwind turbulence intensity and average non-uniformity index at the location of the test section are both less than 2% when the wind speed reaches 6 m/s. Figure 6 illustrates the three-dimensional diagram of TJ-7 wind tunnel. The cross section of the prototype of the bridge deck is shown in the Figure 7. The scale ratio for the sectional model is 1:200. The aerostatic forces of the model are measured at the wind speed of 4–12 m/s to calculate the force coefficients, as shown in Figure 8. Meanwhile, same CFD simulation method as used in Section 4.1 is employed to obtain the force coefficients for comparison purpose.

As shown in Figure 9 (a), the static force coefficients are obtained in the wind tunnel tests in a range of Re = [0.68-2.04]·10⁵ considering three angles of attack [−2°, 0°, 2°]. The sectional model shows obvious Reynolds-number dependency in the drag and lift coefficients. Figure 9 (b) shows that the force coefficients have a good agreement between CFD and wind tunnel test results at $\alpha$ = 0°. Meanwhile, the slopes of the lift and moment coefficients present a satisfactory agreement at $\alpha$ = −2° and 2°. It should be noted that the focus should be on the slopes of the lift and moment coefficients since these values have more decisive influence on the flutter derivatives, as described in Eqs. (8)-(9). Similar comparison can also be found in Montoya et al. (2018b).

Figure 9. Static force coefficients obtained by wind tunnel tests and CFD.

4.3 Cross section property and referenced bridge structure

The change of the closed-box girder shape inevitably affects the natural frequency and modal shape of the bridge structure, which also have significant effects on the flutter vertical wind speed. A real suspension bridge with the span arrangement of 500 m + 1666 m + 500 m utilized by Fang et al. (2022) is introduced in this study. The natural frequency and modal shape of the bridge are then calculated using the finite element model by changing the shape of girder.

As the main components of the top and bottom plate of the steel box girder, U-rib can not only enhance the overall bearing capacity of the steel box girder, but also improve the local bending stiffness and compression stability bearing capacity of the top and bottom plate (Chen & Yang, 2002). For the height and thickness of U-rib, Xiao (2011) suggested the thickness of U-rib should be selected from 6 mm to 10 mm, and the height of U-rib should be selected from 260–320 mm. The equal spacing of the U-rib is adopted in this study as recommended by Yang (2010).

Accordingly, the dimensions of the steel box girder section components are set as follows: The thickness of top plate is 18 mm. The thickness of the bottom plate and inclined webs is 10 mm. The first type of U-rib has the thickness of 8 mm, height of 280 mm and width of 300 mm, which is arranged at the bottom of top plate with equal interval of 300 mm. The second type of U-rib is 6 mm in thickness 260 mm in height and 500 mm in width, which is arranged at bottom plate and bottom inclined webs with equal intervals of 500 and 400 mm, respectively. The I-shape rib has the thickness of 10 mm and height of 120 mm which is arranged at upper and lower inclined webs with equal spacing of 200 mm. Specifically, the cross sections of several cases for the closed-box girders with U-ribs are shown in Figure 10.

Figure 10. Cross sections with ribs.

The details of the finite element model of the bridge are shown in Figure 11. The BEAM4 3D beam element is used to model the stiffening girder, main tower and rigid arms connecting the stiffening girder and suspenders. The main cables and suspenders are simulated using the LINK10 3D rod element. The mass of the girder and second-stage dead load are applied by the MASS21 element.

Figure 11. Finite-element model with a main span of 1666m.

5 Result analysis

5.1 The force coefficients and slopes

The Kriging surrogate model is used to fit the force coefficients and their first-order slopes with respect to the angle of attack of the closed-box girder (Note: the results shown here are processed by the infill sampling criteria). There are 4 parameters (factors) with each of 29 or 43 levels considered in this study. For the convenience of displaying, 25 contour maps are illustrated with interval of 6 levels, as shown in Figure 12.

Figure 12. Force coefficients and slopes.

It can be seen that the drag coefficient mainly depends on the horizontal position (B₁) of the wind fairing point and is less sensitive to the vertical position (H). And it shows an increase trend with the increase of the width of bottom plate (B₂). The lift coefficient increases with B₂. But the sensitivity of B₂ seems more significant that other two parameters. The lift coefficient shows a symmetric pattern with H. However, the axis of symmetry shifts down with the decrease of width of bottom plate (B₂).

As for the moment coefficient, the horizontal and vertical positions of the wind fairing (B₁ and H) have different degrees of importance to it. When the depth of section (D) is small, the lifting moment coefficient almost only depends on the H. The importance of B₁ gradually increases to the same as H with the increase of D. In the other dimension, the moment coefficient shows an increase trend with the increase of the B₂ and D.

The slope of drag coefficient presents a valley domain on the sub-diagonal of each sub-graph with the lowest at upper right corner, namely, the large values of B₁ and H. The peak values are almost located at both ends of the main diagonal when the depth of the section (D) is large. Comparatively, the slopes of lift coefficient and moment coefficient show no obvious variation trend with the shape parameters.

5.2 Flutter derivatives

Based on the results of the above force coefficients and their slopes, the flutter derivatives are calculated using the quasi-steady formulation. To illustrate the variation trend of flutter derivatives with the deck shape parameters, four important flutter derivatives (Matsumoto et al., 1996; Bartoli & Mannini, 2008) $(A_1^{\ast },A_2^{\ast },A_3^{\ast },H_3^{\ast })$ at U* = 6 are shown in Figure 13.

Figure 13. Flutter derivatives (U* = 6).

As can be seen, $A_1^{\ast }$ has a lower value domain (H = 1.45 m-2.05 m and B₁= 50.2 m-55.6 m) when the width of the bottom plate B₂= 27.3 m and the depth of section D = 5.04 m. In addition, $A_1^{\ast }$ has a larger value domain (H = 1.75 m-2.35 m and B₁= 41.8 m-50.2 m) with B₂= 27.3 m and D = 3,72 m. Moreover, $A_1^{\ast }$ changes dramatically and has an obvious left-right or upper–lower symmetry at the design domain. Except for the region near D = 5.04 m, where $A_2^{\ast }$ changes sharply, the changes in other regions are gentle. The distribution of $H_3^{\ast }$ is similar to that of $A_1^{\ast }$ and $A_3^{\ast }$, but its maximum and minimum regions are reversed. The distribution of $A_3^{\ast }$ is almost identical to $A_1^{\ast }$, except that the values of $A_3^{\ast }$ are slightly smaller than $A_1^{\ast }$.

5.3 Dynamic characteristics

Based on the finite element model in Section 4.3, the fundamental frequencies and mode shapes of the bridge can be solved by the Lanczos method. In order to take the modal participation of components such as bridge towers, main cables and suspenders into account during the bridge girder flutter analysis, the equivalent mass and mass inertia moment per meter of the girder is utilized and can be expressed as: (26) $\begin{array}{rl}{m}_{eq}& =\frac{{\int }_{S}m\left(s\right){\varphi }_i^2\left(s\right)ds}{{\int }_G{\varphi }_{h,i}^2\left(x\right)dx}\end{array}$(26) (27) $\begin{array}{rl}{I}_{eq}& =\frac{{\int }_{S}I\left(s\right){\varphi }_i^2\left(s\right)ds}{{\int }_G{\varphi }_{\alpha ,i}^2\left(x\right)dx}\end{array}$(27) where ${\int }_{S}m\left(s\right){\varphi }_i^2\left(s\right)ds$ and ${\int }_{S}I\left(s\right){\varphi }_i^2\left(s\right)ds$ are the modal mass and mass inertia moments of the full bridge for the ith-order vibration mode, which can be normalized as 1. ${\int }_G{\varphi }_{h,i}^2\left(x\right)dx$ and ${\int }_G{\varphi }_{\alpha ,i}^2\left(x\right)dx$ are the integrals of the squares of the ith-order vertical bending and torsional modes of the girder, respectively. Based on the cross-section property of the girder and the bridge structure described in Section 4.3, the frequencies and equivalent mass or mass moment of inertia of fundamental symmetric and antisymmetric vertical bending and torsional modes are calculated for all sampled sections.

For the illustration purpose, the Kriging surrogate model is applied to the frequencies of fundamental symmetric vertical bending and torsional modes, as shown in Figure 14 (a)-(b). The frequency ratio between these two modes, which is considered to have the significant effects on flutter performance (Wang et al., 2014; Abbas & Morgenthal, 2016), is also given in Figure 14 (c). As can be seen, the vertical bending frequency f_h appears to increase with increasing D and decreasing B₁ in Figure 14 (a). Especially, f_h increases as B₂ decreases when D > 5.40 m. And f_h is almost independent of H. Meanwhile, distributions of torsional frequency f_α and frequency ratio $\xi$ is almost same as f_h. The difference, however, is that the value of f_h increases from 0.102 Hz to 0.112 Hz over the entire design domain by approximately 10%, while f_α increases from 0.21 Hz to 0.30 Hz by about 43% and $\xi$ increases from 2.05 to 2.65 by 29%.

Figure 14. The first-order symmetric frequency.

5.4 Cross-validation of initial surrogate model for flutter critical velocity

The flutter critical velocities in the design domain are calculated by using the flutter derivatives in Section 5.2 and dynamic characteristics of the bridge in Section 5.3. The lower value of the flutter critical velocity calculated from the first-order symmetric modes and anti-symmetric modes is used to develop the velocity contour. The accuracy of the Kriging surrogate model is evaluated by concise and effective leave-one-out cross-validation (Lindman, 2012). Relative error $e_i$, and its mean $\overline{e}$ and standard deviation $\sigma$ are selected as evaluation indicators: (28) $e_i=\frac{V_i^C-V_i^K}{V_i^C}\overline{e}=\frac{1}{n}\sum _{i=1}^n|e_i|\sigma =\sqrt{\frac{{\sum }_{i=1}^n{\left(e-\overline{e}\right)}^2}{n-1}}$(28) where $V_i^C$ and $V_i^K$ are the calculation value and Kriging surrogate model prediction value of flutter critical velocity for ith sample point respectively. The results are shown in Figure 15. The maximum and minimum absolute errors are −25.95% and 0.11%, respectively. The mean error $\overline{e}=$5.12% and standard deviation σ = 6.92%.

Figure 15. Relative error of initial Kriging surrogate model leave-one-out cross-validation.

5.5 Model update

This study mainly focuses on the ability of the Kriging model to capture the peak flutter critical velocity. The criteria in Section 2.5 is used to further capture the maximum flutter critical speed until the maximum critical velocity no longer increases, while ensuring that $|e_i|$ < 5%. New samples are added using the hybrid criteria to update the surrogate model. The new experimental design scheme is shown in Figure 16. On the basis of initial 68 samples, 13 groups of shape schemes with the maximum RMSE (No.69-81) and 37 groups of maximum EI, MP and PI (No.82-118) are added.

Figure 16. Uniform experimental design combination of local adding samples scheme.

The prediction error obtained by the leave-one-out cross validation of the updated Kriging surrogate model is shown in Figure 17. RMSE serves as a supplementary test to improve the global accuracy at the boundary of the design domain. These samples added based on RMSE have insignificant effects on improving the accuracy in the peak U_cr region, which are excluded during the calculation of $\overline{e}$ and σ. The maximum and minimum absolute error are 25.95% and 0.07%, respectively. The mean error $\overline{e}=$3.89% and standard deviation σ = 5.71%. Among, the mean error $\overline{e}=$5.12%, standard deviation σ = 6.92% for samples of No.1-68, and the mean error $\overline{e}=$1.62%, standard deviation σ = 1.95% for No.82-118 samples in the peak area. It suggests that the prediction accuracy of the model to achieve an optimal configuration of the bridge deck with maximum flutter critical velocity increases by EI, MP and PI infill sampling criteria.

Figure 17. Relative error of the updated Kriging surrogate model leave-one-out cross-validation.

The flutter critical velocity distributions predicted by the updated surrogate model with interval of 3 levels are shown in Figure 18. As can be seen, when B₂ < 30.1 m, the value of H associated with the peak of flutter critical velocity increases with the increase of D. The value of B₁ at the peak location increases with the increase of D when D > 5.04 m. It indicates that there is an optimal value for the ratio of the wind fairing height to the section depth ((D-H)/D) and the ratio of the bottom plate width to the bridge width (B₂/B₁). Figure 19 are the contour maps of the maximum value in B₂ ∼ D dimensions and B₁ ∼ H dimensions obtained from in Figure 18, respectively. The range of shape parameters of the peak region are B₂= 26.5 m-30.5 m, D = 5.28 m-6.24 m, B₁= 52.2 m-57.8 m and H = 2 m-2.75 m.

Figure 18. Updated flutter critical velocity.

Figure 19. Maximum flutter critical velocity region.

More specifically, the effects of (D-H)/D, B₂/B₁ and three section angles, i.e. θ₁, θ₂ and θ_T on U_cr are shown in Figure 20 and Figure 21. It can be noted that these shape parameters tend to be stable as the increase of U_cr. When the U_cr > 80 m/s, the values of these two ratios of (D-H)/D and B₂/B₁ vary from 0.52–0.64 and 0.49–0.54, respectively. And the ranges of the three section angles are θ₁= 12.5°−18°, θ₂= 12.5°−16.5° and θ_T= 27°−34°. In addition, the average values of these two ratios of (D-H)/D and B₂/B₁ are 0.5865 and 0.5157 while these three angles are θ₁ = 14.90°, θ₂ = 14.33° and θ_T = 29.23°.

Figure 20. Trend of relative height of wind fairing and width of bottom plate (U_cr> 80 m/s).

Figure 21. Trend of θ₁, θ₂ and θ_T (U_cr> 80 m/s).

5.6 Statistical analysis

The updated surrogate model contains a variety of influence information, which can be visually displayed by analysis of variance (ANOVA). Cohen (2013) summarized various evaluation methods for effect and statistical testing ability of ANOVA, and suggested that ${\eta }^{\text{2}}$ can be used as the indicator of the effect: (29) ${\eta }^2=\frac{S{S}_A}{S{S}_T}$(29) where SS_A refers the sum of squares of deviations between groups, SS_T represents the sum of squares for total. Figure 22 (b) shows the ${\eta }^{\text{2}}$ values of these 4 factors and their interactions used in this study.

Figure 22. Main effect and ratio of ${\eta }^{\text{2}}$.

The main effect of these 4 factors on the flutter critical velocity is shown in Figure 22 (a). The red horizontal dashed line is the overall average velocity $\overline{\mathrm{U}}$_cr. When the curve of main effect is parallel to the horizontal axis, it indicates that there is no main effect, i.e. the response value does not change with the value of the factor. The main effect of D on the ${\overline{U}}_{cr}$ has the largest variation, which increase from 55 m/s to 63 m/s at the design domain. The influence of other three factors on $\overline{\mathrm{U}}$_cr is roughly same. But their effects are weaker than that of the girder depth in the design domain. It is worth noting that the trend of effects of H and B₁ on $\overline{\mathrm{U}}$_cr is nearly identical. The interactions of these four factors and their combinations are shown in Figure 22 (b), where * refers combination of the two factors. It can be found that factor D also has the greatest impact on flutter performance, up to 36%. The effects of the other three factors, the second-order and third-order interactions on the flutter vary from 2% to 6%, which verify the conclusion obtained from Figure 22 (a).

Moreover, the comparison of the effects of aerodynamic parameters and dynamic characteristics on the flutter performance due to the change of the deck shape is conducted. Two relative difference coefficients, ${\zeta }_{dynamic}$ and ${\zeta }_{aero}$ are defined as: (30) $\begin{array}{rl}{\zeta }_{dynamic}& =\left(\frac{U_{cr}}{U_{cr.aero}}-1\right)\times 100\mathrm{\%}\end{array}$(30) (31) $\begin{array}{rl}{\zeta }_{aero}& =\left(\frac{U_{cr}}{U_{cr.dynamic}}-1\right)\times 100\mathrm{\%}\end{array}$(31) where $U_{cr.aero}$ and $U_{cr.dynamic}$ are the flutter critical velocities only considering aerodynamic parameters and dynamic characteristics, respectively. The dynamic characteristic of the bridge and flutter derivatives of the bridge deck are fixed using the parameters in the centre point of the design domain for calculating the $U_{cr.aero}$ and $U_{cr.dynamic}$, respectively. The geometric parameters of central sample are B₂= 30.1 m, D = 4.56 m, B₁= 50.2 m, H = 2.1 m. The results are shown in Figure 23.

Figure 23. ability of dynamic characteristics and aerodynamic geometric shape to improve flutter performance.

It can be seen that ${\zeta }_{aero}$ varies from −35% to 45% (46.88 m/s to 85.20 m/s) while ${\zeta }_{dynamic}$ increases from −22% to 27% (43.42 m/s to 81.92 m/s). For the most important parameter D as discussed before, the improvement of ${\zeta }_{dynamic}$ caused by the increase of D is about 30% while ${\zeta }_{aero}$ varies from 0% to 45%. The improvement of the flutter performance by the dynamic characteristics is mainly due to the improvement of the structural stiffness and vibration frequency from the increase of D and B₂ as seen from Figure 23(a) and Figure 14. The ${\zeta }_{dynamic}$ is almost positively correlated with the deck shape parameters while ${\zeta }_{aero}$ shows strongly nonlinear trend. It can be noted that the optimal shape without considering dynamic characteristics is different from that achieved by considering both aerodynamic and dynamic effects as seen from Figure 23(b) and Figure 18(a). The parameters of the optimal shape without considering dynamic characteristics have two regions, i.e. D = 2.88 m - 3.12 m and D = 5.28 m – 6.24 m. For the region of D = 2.88 m - 3.12 m, the bridge deck is very streamlined with favourable aerodynamic performance, but it is usually difficult to meet the structural static and dynamic requirements. It suggests that the change of dynamic characteristics during the aerodynamic optimization of the bridge girder are nonnegligible, especially for the variation of D and B₂.

5.7 Optimal section

Based on above analyses, the cross-section sizes, dynamic characteristic, motion-induced aerodynamic force and flutter critical velocity of EI, PI and MP in the new round are shown in Table 6 and Table 7.

Table 6. Cross-section characteristic of the optimal section.

No.	B2	B1	H	D	$C_{D,0^{\circ }}$	$C_{L,0^{\circ }}$	$C_{M,0^{\circ }}$	$C_D^{\mathrm{\prime }}$	$C_L^{\mathrm{\prime }}$	$C_M^{\mathrm{\prime }}$	$V_i^C$	$V_i^K$	$e_i$
MP	28.7	56.2	2.5	6.00	0.0230	−0.0764	0.0057	−0.0350	1.9098	1.2304	87.10	92.12	−4.8%

Table 7. Dynamic characteristic of the optimal section.

No.	1-A-V-f	1-A-V-m	1-A-T-f	1-A-T-Im	$V_i^C$	1-S-V-f	1-S-V-m	1-S-T-f	1-S-T-Im	$V_i^C$
MP	0.0882	41.83	0.2927	17615	132.2	0.1091	30.23	0.2602	9522	87.1

Note: A, S, V and T represent antisymmetric, symmetric, vertical and torsional modes, respectively. f is the frequency in Hz. The units of mass m and mass moment of inertia per unit length are t/m and t·m²/m.

It can be noted that $|e_i|$ of optimal section is less than 5% and the optimal shape of flutter performance is shown in Figure 24. The shape parameters of the optimal section are B₂ = 28.7 m, D = 6 m, B₁ = 56.2 m, H = 2.5 m, θ₁ = 14.29°, θ₂ = 14.31°, θ_T = 28.61°, B₂/B₁ = 0.51, width-depth ratio = 9.37 with ratio of height of wind fairing to depth of section (D) = 0.58. The flutter critical velocity by CFD is 87.1 m/s. And 145% growth of the critical wind speed is achieved by comparing with the section with the minimum flutter critical velocity of 35.5 m/s in the design domain.

Figure 24. Optimal section.

6 Conclusion

Four geometric parameters of the streamlined closed-box girder, namely the horizonal and vertical coordinates of the wind faring corner point (defined by the width of the section B₁ and vertical distance away from the central point of the bridge deck H) and the intersection point between the inclined lower web and bottom panel (defined by the width of the bottom plate B₂ and depth of the bridge deck D), are defined in this paper to find an optimal configuration considering flutter performance. The effects of the deck shape on the flutter derivatives as well as the dynamic characteristics of the bridge are considered. The uniform design, Kriging surrogate model, CFD and finite element method are combined to achieve the optimization process. The main conclusions are as follows:

1. The use of surrogate model and CFD technique allows a significant reduction of the computation workload as compared with the full factor experiment, i.e. three order of reduction in this study. The efficiency and accuracy can be improved significantly by hybrid infill sampling criteria and preforming the leave-one-out cross-validation.

2. The vertical bending frequency f_h appears to increase with increasing D and decreasing B₁. Especially, f_h increases as B₂ decreases when D > 5.40 m. And f_h is almost independent of H. The distribution patterns of torsion frequency f_α and frequency ratio $\xi$ are almost same as f_h. The difference is that the value of f_h increases from 0.102 Hz to 0.112 Hz over the entire design domain by approximately 10%, while f_α increases from 0.21 Hz to 0.30 Hz by about 43% and $\xi$ increases from 2.05–2.65 by 29%.

3. Optimal ranges for the ratios of (D-H)/D and B₂/B₁ as well as three section angles, i.e. θ₁, θ₂ and θ_T are found to achieve a relatively high flutter critical wind speed (U_cr > 80 m/s). The values of these two ratios of (D-H)/D) and (B₂/B₁) vary from 0.52–0.64 and 0.49–0.54, respectively, with the averages of 0.5865 and 0.5157. The ranges of three angles are θ₁ = 12.5°−18°, θ₂ = 12.5°−16.5° and θ_T = 27°−34° with the mean values of θ₁ = 14.90°, θ₂ = 14.33° and θ_T = 29.23°, respectively. The recommended ranges of these dimensionless parameters provide various favourable shape selections during the preliminary design of the bridge.

4. The depth of the section D is the most significant influence factor on flutter performance accounting for 36%. The effects of other three factors, second-order and third-order interactions all less than 6%. Based on the central sample in the design domain, the increase of D results in that the improvement ratio of the flutter critical wind speed due to the only variation of aerodynamic parameters varies from 0% to 45% nonlinearly. But it also contributes to the enhancement of the flutter performance by 30% when the variation of the dynamic characteristics is only considered. The change of dynamic characteristics in aerodynamic optimization of closed-box girder cannot be ignored, especially for the variations of the depth of section and width of bottom plate.

5. Peak regions of flutter critical velocity are found for the closed-box girder used in this study. That is, optimal sections are achieved. The parameters of the optimal section are B₂ = 28.7 m, D = 6 m, B₁ = 56.2 m, H = 2.5 m, θ₁ = 14.29°, θ₂ = 14.31°, θ_T = 28.61°, B₂/B₁ = 0.51, width-depth ratio = 9.37 with ratio of height of wind fairing to depth of section (D) = 0.58. The flutter critical velocity is 87.1 m/s, which are 145% higher than that of the section with the minimum critical speed.

Acknowledgments

The authors gratefully acknowledge the support of the National Natural Science Foundation of China (52108469, 51978527, 52078383, 52278520), the Fundamental Research Funds for the Central Universities (22120220577).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China [grant numbers 52108469, 51978527, 52078383, 52278520]; Fundamental Research Funds for the Central Universities [grant number 22120220577].