A New Approach for Free Vibration Analysis of Thin-Walled Box Girder Considering Shear Lag Effect

*e thin-walled box girder (T-WBG) is widely applied in the long-span bridge structures during the past decades due to its lighter self-weight and better mechanical properties.*e shear lag effect (SLE), an essential aspect of T-WBGwhich governs the stress and the deformation, is rather necessary to be revealed properly. *e extraordinary issue of T-WBG analysis nowadays is the SLE impact on its dynamical response to external load. *is paper proposes an improved finite element method (FEM) to obtain the realistic vibration characteristics of the T-WBG considering the SLE by theory analysis and formula derivation. Firstly, based on the classical plate and shell theory as well as beam theory, the T-WBGwas divided into shell subunit for the roof and beam subunit for web and floor, respectively. Secondly, a 3-order polynomial which is consistent with the experiment results was adopted as the axial-displacement interpolation function of the roof subunit, whose nodal displacements parameters were also taken as the basic. *irdly, the nodal displacement parameters of the web subunit and floor subunit were deduced by the basic according to the principle of deflection consistency. It is shown through a numerical example that the proposed method is much more economical to achieve reasonable accuracy than traditional FEM analysis software when dealing with the free vibration problem of the T-WBG considering the SLE. Besides, it is also observed that the natural frequency values considering the SLE have a trend of decreasing markedly in general, and the influence of SLE on higher-order frequency is more significant than on the lower one under the boundary condition of cantilever supported, while a contrary effect under the boundary condition of simple supported.


Introduction
SLE was first noticed in aerospace structural engineering practice, which has a direct impact on the detailed structural design of T-WBG accessory structures such as wing spar and flap [1][2][3].In order to reduce self-weight and increase spanning capacity, T-WBG structures nowadays are extensively used in long-span bridges, and the SLE investigation has become extremely important in design and construction of such bridges [4][5][6][7][8][9][10]. is is because (1) the T-WBG used in long-span bridges becomes more flexible as the span becomes longer, leading to a remarkable increase in SLE; (2) the widthspan ratios of the long-span bridges are generally larger nowadays, and theoretical study and experimental observations suggest that the SLE is prone to occur in T-WBG with a large width-span ratio.At present, the T-WBG is widely used in long-span bridges, and its dynamical response to pedestrian load, vehicle load, and so on is much more complex than the static one, which indicates that the vibration analysis of the T-WBG considering the SLE has become a hot spot in the research field.During the past decades, much effort has been directed towards SLE analysis of the T-WBG.Among these methods [11][12][13][14][15][16][17][18][19][20][21][22][23][24], the energy variation method (EVM) and the FEM are considered to be two of the most widely used methods.
For the EVM, numerous studies have contributed to its development.Reissner [25] firstly proposed the EVM for SLE.In his study, a second-order polynomial was used as an approximate axial-displacement interpolation function, trying to describe the actual longitudinal displacement mode considering the SLE.An improved sequential Reissner approach validated by model test and numerical analysis was presented by Guo et al. [26], which used a third-order polynomial as the distribution function.Further research was carried out by Chang [27], in which superposition principle is used to account for the prestress influence on SLE.While these methods above are adopted to deal with SLE problems, difficulties are encountered.ese difficulties are (1) the EVM method is not applicable to section-altered box girder due to the inability for a closed solution achieving; (2) the result reliability is difficult to evaluate and a further revision to the result is not easy; and (3) the EVM theory is complicated such that it is not convenient for an engineering application.
As for the FEM, several methods have been established as an attempt to overcome the problems in the EVM.Jing et al. [28] proposed a FEM for SLE analysis in section-altered box girder with long overhanging flange and demonstrated the capacity of the FEM application to the SLE analysis of section-altered box girder.Tenchev [29] compared the performance of the FEM in SLE analysis with that of model test results and confirmed the results reliability of the FEM.Luo [30] originally presented a FEM-based method called finite segment method, in which a semianalytical finite segment model of planar beam element was established, reducing the three-dimension problem to one-dimension problem, leading to an easier FEM application to SLE analysis.Although some progress has been made in the FEM for SLE analysis of T-WBG, there are mainly two issues remaining.First, the axial-displacement interpolation function, which is directly correlated with the accuracy of the results, is not adequately and reasonably chosen in the FEM.Second, the result accuracy of the FEM is highly dependent on the overall numbers of degrees of freedom (DOFs).Generally speaking, to achieve a valid FEM result, the total element DOF numbers must be "huge and plenty," which leads to a long computation time in modeling and solving.
ird, in order to obtain a relatively accurate solution, an enough dense meshing is customarily conducted, resulting in a great increase of DOFs, which may bring out singularity of the result, sometimes even with wrong conclusions.
To the best of the authors' knowledge, much more studies aiming at solving the static response of the T-WBG considering the SLE have been conducted, while there have been relatively few research studies for the SLE impact on dynamic characteristic of the T-WBG.So, this paper presents an improved FEM for free vibration analysis of T-WBG considering the SLE, in which the roof subunit was modeled as plate element, and the web subunit and floor subunit were simulated as beam elements, respectively.e nodal displacement parameters of the subunits were all only determined by the basic parameters of the roof subunit, resulting in a dramatically decrease in overall element DOFs.A 3-order polynomial was used as an interpolation function for axial-displacement interpolation function, which conforms to the actual axial-displacement mode caused by SLE. e element integral stiffness matrix was deduced on the basis of variation principle.To show the accuracy of the proposed method, the first six natural frequency parameters of the same T-WBG as [31] were calculated for comparison.Besides, an engineering example was studied for the accuracy and efficiency verification of the proposed method, and some valuable conclusions for SLE effect on natural frequency considering the boundary conditions are also achieved.

The Proposed Method
e proposed method is a hybrid approach, consisting of the methods for solving planar stress problem and thin plate bending problem with small deflection. is method is based on three key concepts: (1) the roof subunit is modeled as the plate element whose stress state can be superposed by planar stress state and thin plate bending state with small deflection; (2) the roof subunit is seen as the basic element, whose nodal translation and rotation DOFs are taken as the basic parameters; (3) the web subunit and floor subunit are simulated as beam elements whose displacement parameters can be derived from that of the roof subunit according to the principle of deflections consistency.Figure 1 shows the principle of the proposed method.

e Nodal Displacement Mode of Roof Subunit.
e nodal translation and rotation DOFs of the roof subunit can be written as follows: where ) and u i , v i , ω i , θ xi , θ yi are the nodal axial displacement, transverse displacement, vertical displacement, rotational displacement around x-axis, and rotational displacement around y-axis, respectively, as shown in Figure 2.
In order to analyze the SLE more accurately, a 3-order polynomial conforms to the actual displacement mode was taken as the axial-displacement interpolation function of the roof subunit.For the cross section shown in Figure 3, the axial displacement at any node of the roof subunit can be expressed as follows: where ξ � (2y/k) is the node transverse coordinate in local coordinate system and k is the transverse width of the basic element.
e Lagrange interpolation function is used as the interpolation function of vertical displacement and rotational displacement; the shape function of the element can be expressed by the local coordinate system as follows: Due to space limitations, the elements in the matrix N are not listed here.
e stiffness matrix of the roof subunit in planar stress state can be obtained as follows: where t, d, and k are the thickness, element length in x-axis direction, and element width in y-axis direction, respectively, and B P and D P are the strain matrix and elastic matrix of the roof in planar stress state.

Advances in Civil Engineering
e sti ness matrix of the roof subunit in thin plate bending state with small de ection can be obtained as follows: where B B is the strain matrix of the roof in thin plate bending state with small de ection.

e Nodal Displacement Mode of Web Subunit.
All the plates that constitute the web subunit bear the loads together, and the stress of a single plate is re ected by the beam stress characteristic, which suggests that each plate making up the web subunit can be simulated by beam element.In addition, the web subunit plays a primary part in vertical sti ening on the roof subunit for its high aspect ratio such that only vertical displacement and axial displacement are taken into consideration.
As is shown in Figure 2, each plate making up the web subunit uses a polynomial of the degree one and a cubic polynomial as axial displacement mode and vertical displacement mode, respectively, and the interpolating functions are as follows: where ζ x/d.

Advances in Civil Engineering
Each nodal displacement parameter can be obtained according to the principle of deflections consistency between the roof subunit and web subunit, and the axial displacement parameters of nodes 7 and 8 are obtained on the basis of equations ( 1)-( 3), which are as follows: Similarly, the angular displacements around y-axis of nodes 7 and 8 are as follows: e vertical displacements of nodes 7 and 8 can also be organized by equations ( 1)-( 3), which can be expressed as follows: e left web and right web of the web subunit are all regarded as beam elements, and the angular displacements around y-axis at M end (as shown in Figure 2) are identical with that of nodes 7 and 8.
e centroid angular displacement of the floor subunit around y-axis is linearly interpolated between nodes 5 and 6, which is in line with the two basic assumption of beam element: (1) the longitudinal fiber squeeze of the beam element is not taken into consideration; (2) the rotation angle of each beam section is the same.e angular displacements of the left web and right web around y-axis at M end are given as follows: e axial displacements of the left web and right web at M end are expressed as follows: e axial displacements of nodes 5 and 6 are as follows: e centroid axial displacement of the floor subunit at M end is linearly interpolated between nodes 5 and 6, which is as follows: e vertical displacements of the left web and right web at M end are the same as that of nodes 7 and 8, which can be represented as follows: e centroid vertical displacement of the floor subunit at M end is linearly interpolated between nodes 5 and 6, which is as follows: By using the same method, the displacement mode of each plate at N end (as shown in Figure 2) can also be obtained.
With regard to the left web, let u * l , ω * l be the axial displacement parameter and vertical displacement parameter of the centroid, respectively, which are as follows: e axial displacement and vertical displacement of the left web are as follows: where A is the transformation matrix of u * l and δ, which can be obtained on the basis of equations ( 1), ( 8), ( 13), (16), and (24) and B is the transformation matrix of ω * l and δ, which can also be obtained based on equations (1), ( 10), ( 13), (21), and (25).e elements of matrix A and B are as follows: e remaining elements are zero without exception.As for the right web, let u * r , ω * r be the axial displacement parameter and vertical displacement parameter of the centroid, respectively, which are as follows: e axial displacement and vertical displacement of the right web are as follows: where C is the transformation matrix of u * r and δ, which can be obtained on the basis of equations ( 1), ( 7), ( 14), (17), and (28) and D is the transformation matrix of ω * r and δ, which can also be obtained based on equations (1), ( 9), ( 14), (22), Advances in Civil Engineering and (29).Due to space limitations, the elements of matrix C and D are not listed here.

e Nodal Displacement Mode of Floor Subunit.
To the floor, let u * b , ω * b be the axial displacement parameter and vertical displacement parameter of the centroid, respectively, which are as follows: e axial displacement and vertical displacement of the centroid are as follows: where E is the transformation matrix of u * b and δ, which can be obtained on the basis of equations ( 1), ( 7), ( 8), ( 13), ( 14), ( 18)- (20), and (31) and F is the transformation matrix of ω * b and δ, which also can be obtained based on equations ( 1), ( 11)-( 15), ( 21)-( 23), and (32).Also due to space limitations, the elements of matrix E and F are not listed here.
Based on the displacement modes of the web subunit and floor subunit, the stiffness contributions to element stiffness matrix of the two can be obtained by the minimum potential energy principle, the expression of which is as follows: where  r is the strain energy of the roof subunit and  w and  f are the strain energy of the web subunit and floor subunit, respectively, which can be written as follows: where EA l , EA r , and EA b are the compressive stiffness of the left web, right web, and floor, respectively; EI yl , EI yr , and EI yb are the bending stiffness of the left web, right web, and floor, respectively; K e w and K e f are the stiffness matrix of web subunit and floor subunit, respectively; and F eT is the external load array.
By taking variation on equation (36) and equating it to zero, the stiffness matrix of the web subunit and floor subunit can be obtained as the following: e stiffness matrix of the roof subunit can be obtained by superposition of K e P and K e B , which is as follows: (38)

e Dynamic Equilibrium Equation of FEM.
e vibration system is subjected to three kinds of forces caused by mass, damping, stiffness, and external loads.e dynamic equilibrium equation of FEM can be obtained by dispersing the dynamic system, which is as follows: where M e , C e , and K e are the mass matrix, damping matrix, and the element global stiffness matrix, respectively and € δ, _ δ, δ, and P e are the acceleration vector, velocity vector, displacement vector, and dynamic load vector, respectively.

e Coefficient Matrix.
In actual, the damping has little effect on natural frequency and vibration mode, which indicates that the force induced by damping can be ignored.erefore, only the mass matrix and stiffness matrix are discussed in the following.

e Element Mass Matrix.
e uniform mass matrix is used to analyze the structural vibration, and the method for mass matrix derivation is the same as that of stiffness matrix.Let the material mass density of the element be ρ; the uniform mass matrix of the roof subunit can be obtained according to the shape function from the planar stress state and thin plate bending state with small deflection, which is as follows: where A u is the area of roof subunit.e uniform mass matrix of the web subunit can be obtained in the same way, that is, Advances in Civil Engineering e uniform mass matrix of the floor subunit can be obtained in the same way, that is, where A l , A r , and A f are the area of left web, right web, and floor area, respectively.e mass matrix of the whole element can be obtained by addition of M e r , M e w , and M e f .Due to space limitations, the matrix elements are not listed here.

e Element Stiffness Matrix.
e stiffness matrix of the whole element can also be obtained by addition of K e r , K e w , and K e f , which have been deduced above.Due to space limitations, the matrix elements are not listed here.

2.6.
e Equivalent Nodal Load.If the external load P i (τ) (i � 1, 2, 3, 4) is directly applied to the 4 nodes of the roof subunit, the load vector can be directly substituted into the finite element equation for calculation.In actual, the external load acting on bridge deck often consist of distributed load P(x, y, τ) and concentrated load P(x i , y i , τ); the uniform nodal load array produced by the two kinds of load is as follows: (43)

Numerical Examples and Discussion
In order to demonstrate the accuracy of the proposed method, the first four dimensionless natural frequencies of the same T-WBG as discussed by Mashat [31] were calculated for comparison.e boundary condition of the following analysis is simply supported, and the geometric properties of the roof subunit, web subunit, and floor subunit, as well as material properties, were taken as the same as the reference.e schematic diagram of the square box beam is shown in Figure 4, and the geometry of square box beam is listed in Table 1.e present results are shown in Table 2 together with the reference results.
According to the comparison between the present results and reference results in Table 2, it is shown that the results obtained using the proposed method are close to the reference, which verifies the accuracy of the proposed method.Besides, the frequency results calculated by the proposed method is generally a little greater than the reference ones, which is because in the proposed method, the web subunit and floor subunit are all subjected to tension force, and the rigidity of the two subunits accordingly increases, leading to an increasing of the system stiffness.erefore, the present approach is more accurate.
In addition, a typical engineering example was considered in this section for efficiency and accuracy of the proposed method.
e example is a concrete T-WBG bridge with a main span of 30 m. e cross section layout of the T-WBG is single-box-single-cell with a long cantilever which is 3.55 m long and 0.25 m thick.
e width and thickness of the floor is 7.1 m and 0.25 m, respectively, and the height and thickness of the web is 2 m and 0.4 m.In order to study the effect of different boundary conditions on SLE, two boundary conditions (cantilever supported and simply supported) are adopted in the following research.
e restraint stiffness of the two boundary conditions is fixed by specification parameters of the actual supports.e plan and cross section layouts of the engineering example are shown in Figure 5, and the structural dimensions and material parameters of the engineering example are shown in Table 3.
e SLE analysis of the engineering example was conducted by using the finite element program compiled with the method proposed above (the compiling process of the calculation program is shown in Figure 6).In the program, the bridge was divided into 30 elements in x-direction, and only natural frequencies of undamped free vibration were calculated in order to facilitate the analysis.
e natural frequencies of the preceding 6 modes of the engineering example were compared with that calculated by shell model (shell 181) with ANSYS to verify the efficiency and accuracy of the proposed method.e analysis results are given as Tables 4 and 5.
According to the result comparison between the proposed method and shell model established by ANSYS software in Tables 4 and 5, it is shown that the results obtained using the proposed method are close to that using shell model, which verifies the reliability of the proposed method.Besides, the DOF numbers of the proposed method are reduced by 1922 compared with the ANSYS model, which shows the efficiency of the proposed method.In addition, it is concluded from Tables 4 and 5 that the natural frequency values considering the SLE have a trend of decreasing markedly in general, and the extent of the SLE influence on natural frequency varies with different boundary conditions.For the boundary condition of cantilever supported, the SLE effect on higher-order frequency is more significant than on the lower one, while for the simple supported one, the SLE effect on higher-order frequency is less marked than on the lower one.

Conclusions
Aiming at the structural free vibration characteristics of the T-WBG considering the SLE, an improved FEM has been developed to perform this study.In the proposed method, the roof subunit, which was seen as the basic part of the structure, was analyzed using the methods for solving planar stress problem and thin plate bending problem with small deflection, and the web subunit and floor subunit were all regarded as the attachment structures, the analytical approach of which was the beam element method.e displacement mode for each subunit was derived, and the element global stiffness matrix was deduced on the basis of energy variation principle.e uniform mass matrix and load matrix of the whole element were further derived, and finally the finite element dynamic equilibrium equations were obtained.

Advances in Civil Engineering
A numerical example involving a detailed engineering computational model was presented to demonstrate the e ciency and accuracy of the proposed method.e relevant SLE analysis calculation program was compiled for the computational model above, as well as a nite element model established with ANSYS software for comparison and veri cation of the proposed method.e results show that the proposed method exhibited a better performance in    Note.L/a: the length-to-thickness ratio of the beam; R: reference results; P: present results. 8 Advances in Civil Engineering terms of accuracy and e ciency, as measured by the closeness of the exact values calculated in [31] and by total number of element DOFs compared to the model established with ANSYS.As a result, the proposed method dramatically reduces the number of the DOFs, showing a good ability to improve the calculation e ciency.Moreover, the engineering example demonstrates that the SLE has great in uence on the natural frequency of the T-WBG, and the frequency values with the SLE have a trend of decreasing markedly in general.It is also shown that the extent of SLE e ect on natural frequency varies with di erent boundary conditions.For the boundary condition of cantilever supported, the SLE e ect on higher-order frequency is more signi cant than on the lower one, while for the simple supported one, the SLE e ect on higher-order frequency is less marked than on the lower one.

Figure 1 :
Figure 1: Principle of the proposed method.

Figure 2 :Figure 3 :
Figure 2: Layout of the basic element.

Figure 5 :
Figure 5: Schematic diagram of the engineering example: (a) 3-dimensional graph of T-WBG; (b) cross section layout of T-WBG.

Table 3 :Figure 4 :
Figure 4: Schematic diagram of square box beam: (a) 3-dimensional graph of square box beam; (b) cross section layout of square box beam.

Table 1 :
Geometry of square box beam.

Table 2 :
Comparison of rst four natural frequency parameters in di erent length-to-thickness ratios for square box beam.