Numerical study on surface distributed vortex-induced force on a flat-steel-box girder

ABSTRACT To ensure the safety of bridges, the comfort of pedestrians and vehicles on bridges, the vortex-induced forces on flat-steel-box girders need to be investigated. The total vortex-induced forces integrated on the surfaces of girders have been studied extensively. This study will be mainly focused on the characteristics of surface distributed vortex-induced force, which can be beneficial to the vortex-induced vibration (VIV) reduction with local aerodynamic optimization. Computational Fluid Dynamics (CFD) method is employed for the simulation of the VIVof a flat-steel-box girder. The simulation results are verified through the comparison with the results of corresponding wind tunnel test. A self-adaptive nonlinear fitting method is proposed to determine the vortex-induced force model. A vortex-induced force contribution factor is defined to account the contribution of the components and the surface location. Based on the surface distributed vortex-induced force analysis, several conclusions are made. The vortex-induced force presents a fairly strong multiple-frequency characteristic, and the high-order terms should be carefully considered. Most energy of the VIV comes from the linear and third-order aerodynamic damping terms. In terms of location, the roof takes most of the aerodynamic damping and the total vortex-induced force. The flow separation around the underside upstream corner, the underside downstream corner and the middle downstream comer lead to large third-order aerodynamic damping term. The flow separation around the underside upstream corner has great contribution to the linear aerodynamic stiffness term.


Introduction
Flat-steel-box girders are widely adopted in long-span bridges, for the advantages of light weight, high strength and good aerodynamic performance against flutter instability. Nevertheless, this kind of girder is often subjected to vortex-induced vibration (VIV) at low wind speeds, such as Storebelt Bridge (Larsen, Esdahl, Andersen, & Vejrum, 2000), Trans-Tokyo Bay Crossing Bridge (Fujino & Yoshida, 2002) and Xihoumen Bridge (Laima, Li, Chen, & Li, 2013;Li et al., 2011). Although VIV is a kind of self-limited vibration without disastrous consequences, continual VIVs may cause structural fatigue and discomfort to pedestrians and vehicles. Thus, VIV of flat-steel-box girders should be taken seriously, and an accurate and reasonable prediction is necessary.
It should be noted that most efforts have been focused on the VIVs of circular cylinders, which are easier to reveal some fundamental mechanisms. A large number of studies adopting wind tunnel tests are contained in several reviews, such as Williamson and Govardhan (2004), Sumner (2010), Bearman (2011) and Wu, Ge, and Hong CONTACT Bin Wang wangbinwvb@swjtu.edu.cn (2012). And relative full-scale researches of CFD of circular sections are presented at the early twenty-first century (Guilmineau & Queutey, 2002;Kumar & Dalal, 2006;Tamura, 1999). In fact, compared with circular cylinder, VIV of non-circular section has different characteristics. However, limited research can be found for non-circular sections, especially for the flat-steel-box bridge decks. The experimental results of Diana, Resta, Belloli, and Rocchi (2006) indicated significant nonlinear behaviors of the vortex-induced force on a multi-box deck. Sun, Owen, and Wright (2009) testified that 2D simulation with k-ω RANS turbulence model is appropriate to simulate the wind-induced forces on bluff bodies compared with the 3D simulation. The mechanism and vibration reduction of VIV of a box cross section in the presence of aerodynamic countermeasures were discussed by Sarwar and Ishihara (2010), which showed the influences of fairings and flaps. Hallak et al. (2013) discussed the aerodynamic behavior of a girder in the presence of tall vehicles via 2D CFD model. Section and full aero-elastic model experiments of VIV were conducted and the differences between the two testing results were discussed by Sun, Li, and Liao (2013). However, the mechanism of nonlinear vortex-induced vertical force at resonance stability period on flat-steel-box bridge deck has not been studied extensively.
After lots of exploration about the principles and phenomena of the VIV, only a few researchers proposed mathematical models which are suitable for long-span bridges (Ehsan & Scanlan, 1990;Larsen, 1995;Scanlan, 1981;Simiu & Scanlan,1996). However, these models cannot reflect the features of vortex-induced force comprehensively, and the traditional identification methods introduced in those studies cannot identify precise model parameters to reconstitute the displacement response which could agree well with the CFD results. Aimed at the defects of Scanlan's model and traditional parameters identification method, lots of efforts have been made and some new opinions have been proposed (Meng, 2013;Zhu, Meng, & Guo, 2013). They predicted the responses of flat closed-box cross sections using Scalan's model, and it was found that this model was not appropriate for flat closed-box sections. Afterwards, a new nonlinear model of vortex-induced force for flat closed-box deck, which could reconstitute the vortex-induced force more precisely, was proposed by Zhu et al. (2013). Wind tunnel tests were adopted in the studies mentioned above, which might be influenced by many unpredictable factors, even the results might differ from each other at different wind tunnels with the same research targets (Sarkar, Caracoglia, Haan, Sato, & Murakoshi, 2009). With the development of the computational techniques, numerical simulation gradually becomes a powerful tool to investigate the VIV of bridges, for its convenience and quickness.
In this study, CFD was utilized to explore the vortexinduced force of a flat-steel-box girder. In Section 2, the simulation method and numerical model were introduced firstly. The results obtained from CFD including both the displacement and the vortex-induced force were then compared with those in wind tunnel test in Section 3, to validate the accuracy of the numerical simulation method. In Section 4, a Self-adaptive Nonlinear Fitting Method (SANFM) was introduced, in order to establish an appropriate mathematical model for vortexinduced force of the flat-steel-box girder. The corresponding model was compared with an advanced vortexinduced force model proposed in Zhu et al. (2013). Based on SANFM, the vortex-induced force model of the flat-steel-box girder from the CFD results was then obtained. A contribution factor based on energy spectrum was defined, and the contributions of vortexinduced force components and the locations on the girder to total vortex-induced force were analyzed in Section 5. Finally, some conclusions of this study were drawn in Section 6.

Simulation scheme
The flows around the girder should be simulated to obtain the vortex-induced force. The governing equations of incompressible viscous fluid include mass conservation equation (continuity equation), momentum conservation equation (Navier-Stokes equation) and energy conservation equation. For the air flow around the girder, the heat flux is very limited, thus the energy equation is ignored when calculating VIV responses. After being Reynolds averaged, the governing equations become where ρ and μ are the density and the dynamic viscosity coefficient of air, respectively; subscript i, j represent the coordinate axis in space with i = 1, 2 and j = 1, 2 for two-dimensional simulation; x i is the coordinate value along the i axis; u j and u j are the average velocity and the fluctuating velocity along j axis, respectively; t is the time; p is the pressure; and −ρu i u j is the so-called Reynolds stress. In this study, the commercial software FLUENT was employed with SST k-ω turbulence model to model the Reynolds stress. Second-order scheme was used to discretize the governing equations, based on the finite volume method. SIMPLEC algorithm was employed for the pressure-velocity coupling. The least squares cellbased method was utilized for the gradient term. The pressure was discretized in a standard way. The momentum, turbulent kinetic energy and specific dissipation rate were discretized in second-order form. The above scheme was applied effectively in Wang, Xu, Zhu, Cao, and Li (2013).

Computational domain, boundary conditions and meshing
A flat-steel-box girder of a real long-span bridge is taken as an object, and the wind tunnel test results have been discussed in Zhu et al. (2013). The Wind tunnel tests were carried out in the TJ-3 boundary layer wind tunnel. This wind tunnel is 15 m wide, 2 m high and 14 m long. The turbulence intensity is less than 2.5%. For the sectional model, the scale ratio is 1:20, and the model length is 3.6 m. Moreover, an improved force measurement technique for vortex-induced force (VIF) was adopted in the wind tunnel tests. In order to validate the numerical simulation, the wind attack angle, wind speed and scale ratio in this study are set as the same as these in the wind tunnel tests. That is: the wind attack angle is+5°, the wind speed is 9.1 m/s, and the scale ratio is 1:20. The corresponding Reynolds number is 1.06 × 10 5 based on the height of the girder.
The computational domain is shown in Figure 1. B and D are the width and the depth of the scaled girder, respectively. The width and the depth of the whole computational domain are 17B and 28D, respectively, with a blocking rate of 3.57%. The computational domain is divided into Rigid Domain, Dynamic-Mesh Domain and Stationary Domain as in Huang, Liao, Wang, and Li (2009). Rigid Domain has a rectangle area of 1.12B × 1.25D. The distances from the left edge, right edge, upper edge and lower edge of this domain to the corresponding outermost edges of the girder are 0.05B, 0.1B, 0.125D and 0.125D, respectively. In order to ensure the grid quality near the deck, a body-fitted structured quadrilateral mesh with 15 boundary layers is generated in Rigid Domain. Along the deck, the grid sizes are arranged between 9.6 mm ∼ 18.8 mm, with a majority of 10 mm, which is about 1/3200 of the total width of the girder. Around the railing, the grid sizes are near 10 mm, which is about 1/17.5 of the width of the railing. Around the track maintenance track, the grid sizes are about 4 mm, which is about 1/30 of the width of the maintenance track. The mesh in Dynamic-Mesh Domain is filled with triangular grids. This domain is 3B in width and 10D in height. Stationary Domain is extended from the Dynamic-Mesh Domain to the out boundaries. For the flow around the girder, a vortex street may appear at the downstream region behind the girder in Stationary Domain. The grids in this domain should be densified and identified as Wake Domain to simulate this phenomenon clearly. Wake Domain is a longnarrow quadrilateral of 9.5B × 10D, filled with square grids. Thus, Stationary Domain contains Wake Domain and Exterior Domain. Figure 1 also shows the boundary conditions. The wind attack angle is produced by rotating the inlet flow while the girder keeps horizontal. The left out edge and the lower out edge are the sources of the wind, so they are defined as velocity inlets with a speed of 9.1 m/s. And a turbulent intensity of 0.02 and a turbulent viscosity ratio of 2 are assigned to them. After passing the computational domain, the wind flows out through the right out edge and upper out edge. Thus, the right edge and the upper edge are specified as pressure outlets with a gauge pressure of zero. The whole girder is defined as wall condition. In the whole domain among these boundaries, the flows are solved based on Eqs. (1-2).

Movement simulation of girder
Aerodynamic lift force F t is calculated through CFD and applied on the girder. The movement of the girder in the vertical direction is controlled by a basic dynamic equation as where y t ,ẏ t andÿ t are the vertical displacement, velocity and acceleration of the girder at time t, respectively; m, c and k are the mass, mechanical damping coefficient and stiffness of the girder, respectively. m and k are converted from the parameters in the wind tunnel test of Zhu et al. (2013), which are 50.61 kg/m and 15742.76 N/m, respectively. Moreover, in order to improve the computational efficiency when verifying the independence of grid and time step, the damping coefficient is a constant of 8.926 N·s/m. Fourth-order Runge-Kutta method is used to solve Eq. (3) to obtain the responses of the girder at each time step. The vortex-induced force at the previous time step is applied approximately as F t . From the results in Section 3, this approximation can be accepted. Dynamic mesh approach was utilized in this study according to Batina (1990) and Huang et al. (2009). Rigid Domain (see Figure 1) is linked with the girder and moves with it together. The grids in Dynamic-Mesh Domain are deformed at every time step. And the elastic smoothing method together with local remeshing method is applied to the deformation of the grids in Dynamic-Mesh Domain. When the deformation is very small, the elastic smoothing method is applied. The grid lines are similar to springs, any ending node moves elastically. The elastic rigidity k ij of the two conterminous nodes i and j is where x i and x j are the locations of node i and node j. The whole elastic force imposed on any node should be zero at any time, that is where n i is the number of the nodes near to node i. Thus, in Dynamic-Mesh Domain, the aforesaid equation can be defined as where K is the elastic rigidity of this virtual spring between nodes; and X is the displacement vector of the nodes.
When the deformation of the grids is relatively large, the elastic smoothing method is not applicable, but the local remeshing method, introduced in Fluent 6.3 user's guide (2006), is applied. This method includes four steps. First, all the grids whose cell skewness or sizes exceed the formulary range are marked. Second, all the marked grids are deleted and empty cavities are formed around these grids. Thirdly, these empty cavities are reconstructed with new grids. The quality of the new grids is better than that of the marked grids. Finally, the values of these physical quantities on the new grids are interpolated from the former grids. The principle of the interpolation is to guarantee that the flow physical parameters are consistent before and after interpolation.

Independence study of grid and time step
Three time steps of 0.000223s, 0.000445s and 0.00089s are chosen to check the independence of the time step on the numerical results, corresponding to T/1600, T/800 and T/400, respectively (T is the natural vibration period of the girder). For each case, the results in the first 10 seconds are ignored, and the results in the following 10 seconds are chosen to be compared. Table 1 shows the computed results at different time steps. As can be seen from Table 1, the amplitude of the vertical displacement enlarges with the increase of the time step. And the Strouhal number reduces with the increase of time step. The results of the first two time steps (T/1600 and T/800) are somewhat similar. Compared with other time steps, it can be seen that the time step of 0.000445s is precise   enough, thus, it is adopted in the following numerical simulation.
In the process of the numerical simulation, the y + value has great influence on the calculation results, which is decided by the height of the first layer grids. Consequently, as can be seen in Table 2, three meshing schemes with different heights of the first layer grids are created to check the independence of the numerical results on the grid sizes. G1 has the largest grid number of 0.35 million, the height and maximum y + of its first layer grids are 0.0015 mm and 0.2, respectively. G2 increases the height of the first layer grids to 0.0075 mm, with 0.24 million grids. The maximum y + value is 0.7. Its grid distribution is shown in Figure 2. G3 further increases the height of the first layer grids to 0.015 mm, and the maximum y + value is 1.5, with 0.14 million grids. The computed results on different grids are shown in Table 2. The results of G1 are very close to those of G2, but the results of G3 are far away from those of G1 and G2. In consideration of both the precision of the numerical simulation and the computational efficiency, the meshing scheme G2 is finally applied to the following numerical simulation.

Simulation results and validation
In fact, the structural damping ratio may vary nonlinearly with the change of vibrating amplitude. Zhu et al. (2013) simplified the nonlinear variation of the structural damping ration with a linear equation: where ξ is the structural damping ratio; y is the vibrating amplitude. In this numerical study, the variation of structural damping is also taken into consideration. It is hard to input the structural damping timely in the simulation due to the fact that the vibrating amplitude needs to be determined timely at first. Actually, the moments of the amplitudes corresponds to the moments of the maximum of the potential energy and minimum of the kinetic energy. There is a quadratic dependence of the maximum amplitude on the total energy. Therefore, the linear relationship between the structural damping ratio and the vibrating amplitude is equivalent to a quadratic relationship between the structural damping ratio and the total energy. Then Eq. (7) can be written as where E is the sum of the potential energy and kinetic energy, which can be easily determined timely in the simulation.

Displacement response
As shown in Figure 3(a), the simulated displacement time history of the girder presents the phenomena of self-excited and amplitude limiting, which are the typical characteristics of VIV. Figure 3(b) compares the displacement response of numerical simulation with that of wind tunnel test  at resonance stability period.
In Figure 3(b), for the convenience of comparison, n is determined as a certain moment at the resonance stability period. It can be seen that the displacements of the numerical simulation agree well with the experimental results in amplitude. Figure 3(c) compares the displacement response spectrogram of numerical simulation with that of wind tunnel test  at resonance stability period. The dominant frequency can be captured well in the simulation, but the components of other frequencies are not predicted. Table 3 shows the representative results of the numerical simulation and the wind tunnel test of Zhu et al. (2013). The dominant frequency and the amplitude of the two methods are close to each other with very little differences.

Vortex-induced force
In the wind tunnel tests , the fluctuating force obtained through the balance contains four parts: the pure vortex-induced force on the oscillating girder, the inertia force of the model, the non-wind-induced additional inertia force, and the non-wind-induced additional damping force. The non-wind-induced additional force (including both the inertia and damping components) can be represented as whereF 0 v (t) is the non-wind-induced additional force; y,ÿ are the velocity and the acceleration of the girder; m 0 and c 0 are the non-wind-induced additional mass and damping coefficient, respectively. m 0 and c 0 can be identified though free vibration of the girder with zero wind speed condition.
Comparatively, the fluctuating force in the simulation contains the pure vortex-induced force, the nonwind-induced additional inertia force, and the nonwind-induced additional damping force. Consequently, the vortex-induced force needs to be calculated through subtracting of the non-wind-induced additional inertia and damping force from the total force, as whereF V (t) is the fluctuating part of the lift force; and F VI (t) is the vortex-induced force. With a free vibration in the zero wind speed condition,F V (t) = 0. As a result, m 0 and c 0 can be identified directly from Eq. (10). They are 6.051 kg and 3.098 N·s/m in this study, respectively. Figure 4(a-c) display the time histories ofẏ,ÿ,F V at the resonance stability period in the simulation. The vortex-induced forces of the girder are obtained through Eq. (10) and shown in Figure 4(d). The vortexinduced force and amplitude spectrum, calculated via CFD results, are compared with the experimental results , as shown in Figure 4(e,f). It can be seen that the vortex-induced force from the simulation and the wind tunnel tests agree with each other. The one calculated from the wind tunnel test is not very smooth maybe due to the measurement error, while the one calculated by CFD is smooth and stable. The vortex-induced forces in both wind tunnel test and CFD present a fairly strong multiple-frequency characteristic, which can be interpreted as a symbol of nonlinear behavior of the vortex-induced force.

Self-adaptive nonlinear fitting method
A reasonable and reliable mathematical model for the vortex-induced force is the basic premise for the analysis of vortex-induced vibration. Vortex-induced force on a structure depends on the displacement y, time derivativesẏ of the structure. In Scanlan's empirical linear model (Simiu & Scanlan, 1996), linear terms of y anḋ y are involved. In Scanlan's empirical nonlinear model (Simiu & Scanlan, 1996), a nonlinear term Y 1 ε( y D ) 2ẏ U beside the linear terms of y andẏ is adopted. Zhu et al. (2013) improved Scanlan's empirical nonlinear model by adding another nonlinear termY 3ẏ U y D . Considering all the linear and nonlinear terms, a general form for the vortex-induced force can be expressed as where ρ is the air density; U is the wind speed; D is the height of the girder; ω is the circular frequency of vortex shedding; is the phase difference between the vortex shedding and displacement response; i and j are the order ofẏ and y, respectively; P ij are the coefficients. In order to determine the terms of the vortex-induced force model in Eq. (11) for the flat-steel-box girder appropriately, a Self-adaptive Nonlinear Fitting Method (SANFM) is proposed. In this process, nonlinear least square fitting based on Newton-Raphson method is utilized to get a fitted vortex-induced force model F m VI . The residue between F m VI and the vortex-induced forceF VI can be expressed as where n is the number of data points. The object of the self-adaptive nonlinear fitting process is to get acceptable vortex-induced force model with small residue R. And F VI (t i ) is fitted to get F m VI at each time step. At first, F m VI is obtained without the term i=∞,j=∞ i=0,j=0 P ijẏ i y j in Eq. (11). A first residue R 0 is calculated using Eq. (12). Later on, progressively increasing i and j with increment of 1 in the expression of F m VI in Equation (11). New F m VI is obtained by nonlinear least square fitting, and the residue of the k th pair of (i, j) k is calculated as R k . Also, a criterion is developed to judge whether the k th pair of (i, j) k is acceptable. In this study, the relative residue R k is taken as In this study, the critical value of R k is 1%. If R k > 1%, (i, j) k is accepted and the i th order ofẏ and j th order of y are selected in the vortex-induced force model. If R k < 1%, (i, j) k is not accepted. Through the above process, the vortex-induced force model can be obtained in a self-adaptive way.

Effectiveness of SANFM
The typical vortex-induced vortex models for bridges, e.g. Scanlan (1981)and Larsen et al. (2000), were proposed mostly based on the responses of vortexinduced vibration, aiming at reflecting the response characteristics. Thus, they are closing processes without the participation of vortex-induced forces, which may result in inaccuracy prediction of the vortex-induced forces and the displacement responses. To fulfill the ends, an advanced nonlinear mathematical model was proposed by Zhu et al. (2013)as follows   Table 4. The final vortex-induced force model is Through SANFM, more 2nd and 3th terms appear compared with the model in Eq. (14). Figure 5(a) shows the comparison between the vortex force predicted using Eq. (15) and Eq. (14) and the tested one. It can be seen that Zhu's model cannot predict the vortex-induced force accurately at the position of peaks, and the model determined through SANFM is more precise to reconstitute the vortex-induced force.

Vortex-induced force model based on CFD
As in Section 4.2, the effectiveness of SANFM is validated based on the tested vortex-induced force. In this section, the vortex-induced force model based on the simulation is determined. The process of SANFM is shown in Table 5. The final mathematical model of vortex-induced force for the flat-steel-box girder is Compared with Eq. (15), P 13ẏ U y 3 D 3 is added. From Table 5, the vortex-induced force model using Eq. (15) has a very low relative residue (R k /max. F VI ) of 0.02. The additional term P 13ẏ U y 3 D 3 decreases the residue nearly to zero. Table 6 shows the model parameters in Eq. (16). Figure 6(a) displays the comparison of the simulated and reconstituted vortex-induced vertical force of this model. It can be seen that the reconstituted vortexinduced force is in conformity with the simulated one. Figure 6(b,c) show the comparison between the simulated displacement through CFD and the solved displacement response using the vortex-induced force model. It can then be found that these two time histories agree well with each other. It can be concluded that P 10 , P 01 , P 20 , P 11 , P 02 , P 30 , P 21 , P 03 0.278 1.535 0.020 YES 10 P 10 , P 01 , P 20 , P 11 , P 02 , P 30 , P 21 , P 03 , P 40 9.380 −3274.161 0.688 NO 11 P 10 , P 01 , P 20 , P 11 , P 02 , P 30 , P 21 , P 03 , P 31 8.485 −2951.978 0.622 NO 12 P 10 , P 01 , P 20 , P 11 , P 02 , P 30 , P 21 , P 03 , P 22 0.838 −201.499 0.061 NO 13 P 10 , P 01 , P 20 , P 11 , P 02 , P 30 , P 21 , P 03 , P 13 0.100 64.089 0.007 YES the above results demonstrate the reliability and accuracy of the new mathematical model of vortex-induced force through SANFM for the flat-steel-box girder. In this model: the components P 10ẏ D and P 10 y D are the linear aerodynamic damping force and the linear aerodynamic stiffness force respectively; the pure high-order speed components P 20ẏ 2 U 2 , P 30ẏ 3 U 3 represent the nonlinear aerodynamic damping forces; the pure high-order displacement components P 02

Contribution factor
The components of the vortex-induced force are decomposed from each term in the expression of the model in Eq. (16). In order to analyze the effects of each component, a contribution factor is defined based on energy  spectrum. Fourier transform is firstly used for the signal processing, as shown in Eq. (17).Thus, the amplitude spectrums of each component and the total vortexinduced force are obtained. Then the energy spectrums are obtained, which are the square of the amplitude spectrums. The contribution factor is defined as C O , as shown in Eq. (18).
where j is the number of each component and ranges from 1 to 10; F aj (t) is the j th component of the vortexinduced force; f aj(ω) and f V(ω) are the spectrum of the j th component and the total vortex-induced force, respectively; and C Oj (ω) is the contribution factor of the j th component at frequencyω.

Contribution of components
Through the contribution factor introduced above, the contributions of the ten components in Eq. (16) to three main frequencies (Base-frequency, Double-frequency, and Treble-frequency) are obtained, as shown in Table 7. It can be seen that the dominant frequency of the firstorder terms (P 10 , P 01 ) is 2.718 Hz, equal to the natural frequency of the girder. The dominant frequency of the second-order terms (P 20 , P 11 , P 02 ), the third-order terms (P 30 , P 21 , P 03 ), and the fourth-order term (P 13 ) is the double, the triple, and the four times of the natural frequency. The vortex-shedding force (V S sin(ωt + )) corresponds to the natural frequency. It also can be found that the second-order terms (P 20 , P 11 , P 02 ) take a considerable proportion of 93.3% to the double-frequency of 5.422 Hz. The third-order terms (P 30 , P 21 , P 03 ) take a significant proportion of 99.6% to the triple-frequency and also take a remarkable proportion of 18.0% to the basefrequency of 2.713 Hz. It means that the high-order frequencies can be presented by high-order terms, and the responses of the vortex-induced force with high-order frequencies can be reconstituted via the model with highorder terms. Furthermore, the vortex-shedding force (V S sin(ωt + )) hardly contributes to high-order frequencies, but takes a proportion of 5.2% to the basefrequency. It can be concluded that the vortex-shedding frequency approximately equal to the base-frequency. Force-displacement phase plane is adopted (Figure 7) to analyze the energy evolution of each vortex-induced force component at one cycle of resonance stability period. The vertical displacement is defined as x-axis, the vortex-induced force is taken as y-axis.
From Figure 7, the linear aerodynamic damping term P 10ẏ U forms an ellipse along the clockwise direction, so it inputs energy to VIV. The linear aerodynamic stiffness term P 01 y D forms a straight line across the second and the fourth quadrants, doing no work to VIV. The second-order aerodynamic damping term P 20ẏ 2 U 2 forms a vertical-axis symmetric parabolic at the third and the fourth quadrants, so it has nothing to do with the energy. The second-order cross term P 11ẏ U y D forms a biaxial symmetry figure-of-eight curve, and it inputs energy to VIV with positive displacement and dissipate energy when the displacement is negative, but the energy keeps balance within a cycle.
The second-order aerodynamic stiffness term P 02 y 2 D 2 forms a vertical-axis symmetric parabolic at the first and the second quadrants, it also has nothing to do with the energy. The third-order aerodynamic damping term P 30ẏ 3 U 3 forms a closed curve like a nucleus of jujube, and it inputs energy to VIV. The other two third-order terms P 21ẏ 2 U 2 y D and P 03 y 3 D 3 are both centrosymmetric curves, doing no work. The fourth-order cross term P 13ẏ U y 3 D 3 forms a similar curve with the second-order cross term P 11ẏ U y D , so it has the same characteristics with P 11ẏ U y D . The vortex-shedding force (V S sin(ωt + )) forms a centrosymmetric ellipse across the second and the fourth quadrants with quite small values, so it inputs little energy to VIV. It can be found that, at the resonance period of VIV, the linear and the third-order aerodynamic damping forces input most of the energy to VIV.

Contribution of location
The girder without auxiliary takes 99.2% in terms of the contribution to vortex-induced force. Thus, this study   focuses on the surface of the girder without auxiliary elements. In order to analyze the contribution of different location, the surface of the girder is divided into eight parts (see Figure 8). They are the underside of upstream web (UnUW), the upstream floor (UF), the downstream floor (DF), the underside of downstream web (UnDW), the upper of downstream web (UpDW), the downstream roof (DR), the upstream roof (UR) and the upper of upstream web (UpUW). Table 8 shows the contribution factors of each part to the vortex-induced force components and the total vortex-induced force. Figure 9 shows the streamline and the Q-isolines (Hunt, Wray, & Moin, 1988) around the girder in a time period T. The roof (including both upstream and downstream roofs) makes up 79.03% of the total vortex-induced force, which can also be reflected in Figure 9. In Figure 9, a main vortex forms on the upstream roof, moves along the roof to the downstream roof, and sheds from the end of the downstream roof. The vortex development on the roof is the main source of the vortex-induced force of the flat-steel-box girder. The vortex-induced force on the underside of downstream web takes 9.17% of the total force, which is caused by another vortex developed under the downstream web (see Figure 9). As described in Section 5.2, the linear and the third-order aerodynamic damping forces P 10ẏ U and P 30ẏ 3 U 3 contribute to the total aerodynamic damping. Like the total vortex-induced force, the roof and the underside of downstream web take most of the aerodynamic damping. Meanwhile, the downstream web makes more than 5% contribution to the aerodynamic damping. The linear, third-order, and the fourth-order aerodynamic stiffness forces P 01 As the flow might separates at the surface corners of the girder, another surface division including corners is shown in Figure 10. The upper, middle, and underside upstream corners are arranged in Parts UpUC, UC, UnUc, respectively. The upper, middle, and underside downstream corners are set in Parts UpDC, DC, and UnDC, respectively. The middle parts of the roof and floor without any corner are chosen as Part MR and MF, respectively. Table 9 shows the contribution factors to the vortex-induced force and components with this surface division. The vortexes on the roof mainly are generated from the upper upstream corner and influenced by the upper downstream corner. In this regard, Parts UpUC, MR, and UpDC are counted together and take 80.45% of the vortex-induced force, which is consistent with the observation that the roof makes up the most vortex-induced force. From the flow development in Figure 9, the underside upstream corner, the underside downstream corner, and the middle downstream corner take responsibility for the flow separation. It can be seen in Table 9 that Part UnUC and Part MF (corresponding to the underside upstream corner) provides 27.04% of the linear aerodynamic stiffness term and 12.58% of the third-order aerodynamic damping term. Part UnDC (corresponding to the underside downstream corner) takes 13.64% of the third-order aerodynamic damping term and more than 5% of the linear aerodynamic damping term and the high-order aerodynamic stiffness term. Part DC (corresponding to the middle downstream corner) makes up the 15.08% of the third-order aerodynamic damping term and more than 5% of the high-order aerodynamic stiffness term.

Conclusions
Studies on the vortex-induced force of a flat-steel-box girder are carried out through CFD simulation with SST k-ω turbulence model in this study. A CFD numerical model is set up. Using the local remeshing method and elastic smoothing method, the movement of the girder is achieved. Through the independence studies, an effective time step and mesh scheme is determined. The simulation results are validated with the previous wind tunnel tests, it could be found that the numerical simulation method is reliable with high accuracy. The vortex-induced force on the flat-steel-box girder presents a fairly strong multiple-frequency characteristic.
A Self-adaptive Nonlinear Fitting Method (SANFM) is proposed based on Newton-Raphson method. The effectiveness of SANFM is validated on the tested vortexinduced force in the wind tunnel. Then the mathematical model of the vortex-induced force for the flat-steel box girder based on CFD is determined via SANFM. To account the effects of the components in the model and the effects of surface location to the vortex-induced force, a contribution factor is defined based on the energy spectrum. The linear and the third-order aerodynamic damping terms input most of the energy to the VIV. The roof of the flat-steel-box girder makes up 79.03% of the total vortex-induced force, and the vortex development on the roof is the main source of the vortexinduced force. The roof takes most of the aerodynamic damping terms including the first linear and the thirdorder terms. The flow separation from the underside upstream corner, the underside downstream corner and the middle downstream comer has large contribution to the third-order aerodynamic damping. Besides the roof, the flow separation from the underside upstream corner has great contribution to the linear aerodynamic stiffness term, while the flow separation near the underside downstream corner and the middle downstream corner has certain contribution to the high-order aerodynamic stiffness term.

Disclosure statement
No potential conflict of interest was reported by the authors.