#### Abstract

The 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. The shear lag effect (SLE), an essential aspect of T-WBG which governs the stress and the deformation, is rather necessary to be revealed properly. The extraordinary issue of T-WBG analysis nowadays is the SLE impact on its dynamical response to external load. This 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-WBG was 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. Thirdly, 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.

#### 1. 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–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–10]. This 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 width-span 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–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. These 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. Third, 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. The 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. The 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.

#### 2. The Proposed Method

The proposed method is a hybrid approach, consisting of the methods for solving planar stress problem and thin plate bending problem with small deflection. This 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.

##### 2.1. The Nodal Displacement Mode of Roof Subunit

The nodal translation and rotation DOFs of the roof subunit can be written as follows:where and , , , , 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 is the node transverse coordinate in local coordinate system and *k* is the transverse width of the basic element.

The 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.

The 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 and are the strain matrix and elastic matrix of the roof in planar stress state.

The stiffness matrix of the roof subunit in thin plate bending state with small deflection can be obtained as follows:where is the strain matrix of the roof in thin plate bending state with small deflection.

##### 2.2. The 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 reflected 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 stiffening 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 .

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:

The vertical displacements of nodes 7 and 8 can also be organized by equations (1)–(3), which can be expressed as follows:

The 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. The 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. The angular displacements of the left web and right web around *y*-axis at M end are given as follows:

The axial displacements of the left web and right web at M end are expressed as follows:

The axial displacements of nodes 5 and 6 are as follows:

The centroid axial displacement of the floor subunit at M end is linearly interpolated between nodes 5 and 6, which is as follows:

The 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:

The 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 , be the axial displacement parameter and vertical displacement parameter of the centroid, respectively, which are as follows:

The axial displacement and vertical displacement of the left web are as follows:where **A** is the transformation matrix of and , which can be obtained on the basis of equations (1), (8), (13), (16), and (24) and **B** is the transformation matrix of and , which can also be obtained based on equations (1), (10), (13), (21), and (25). The elements of matrix **A** and **B** are as follows:

The remaining elements are zero without exception.

As for the right web, let , be the axial displacement parameter and vertical displacement parameter of the centroid, respectively, which are as follows:

The axial displacement and vertical displacement of the right web are as follows:where **C** is the transformation matrix of and , which can be obtained on the basis of equations (1), (7), (14), (17), and (28) and **D** is the transformation matrix of and , which can also be obtained based on equations (1), (9), (14), (22), and (29). Due to space limitations, the elements of matrix **C** and **D** are not listed here.

##### 2.3. The Nodal Displacement Mode of Floor Subunit

To the floor, let , be the axial displacement parameter and vertical displacement parameter of the centroid, respectively, which are as follows:

The axial displacement and vertical displacement of the centroid are as follows:where **E** is the transformation matrix of 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 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 is the strain energy of the roof subunit and and are the strain energy of the web subunit and floor subunit, respectively, which can be written as follows:where , , and are the compressive stiffness of the left web, right web, and floor, respectively; , , and are the bending stiffness of the left web, right web, and floor, respectively; and are the stiffness matrix of web subunit and floor subunit, respectively; and 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:

The stiffness matrix of the roof subunit can be obtained by superposition of and , which is as follows:

##### 2.4. The Dynamic Equilibrium Equation of FEM

The vibration system is subjected to three kinds of forces caused by mass, damping, stiffness, and external loads. The dynamic equilibrium equation of FEM can be obtained by dispersing the dynamic system, which is as follows:where , , and are the mass matrix, damping matrix, and the element global stiffness matrix, respectively and , , , and are the acceleration vector, velocity vector, displacement vector, and dynamic load vector, respectively.

##### 2.5. The 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. Therefore, only the mass matrix and stiffness matrix are discussed in the following.

###### 2.5.1. The Element Mass Matrix

The 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 is the area of roof subunit.

The uniform mass matrix of the web subunit can be obtained in the same way, that is,

The uniform mass matrix of the floor subunit can be obtained in the same way, that is,where , , and are the area of left web, right web, and floor area, respectively.

The mass matrix of the whole element can be obtained by addition of , , and . Due to space limitations, the matrix elements are not listed here.

###### 2.5.2. The Element Stiffness Matrix

The stiffness matrix of the whole element can also be obtained by addition of , , and , which have been deduced above. Due to space limitations, the matrix elements are not listed here.

##### 2.6. The Equivalent Nodal Load

If the external load 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 and concentrated load ; the uniform nodal load array produced by the two kinds of load is as follows:

#### 3. 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. The 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. The schematic diagram of the square box beam is shown in Figure 4, and the geometry of square box beam is listed in Table 1. The present results are shown in Table 2 together with the reference results.

**(a)**

**(b)**

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. Therefore, 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. The example is a concrete T-WBG bridge with a main span of 30 m. The 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. The 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. The restraint stiffness of the two boundary conditions is fixed by specification parameters of the actual supports. The 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.

**(a)**

**(b)**

The 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. The 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. The 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.

#### 4. 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. The displacement mode for each subunit was derived, and the element global stiffness matrix was deduced on the basis of energy variation principle. The uniform mass matrix and load matrix of the whole element were further derived, and finally the finite element dynamic equilibrium equations were obtained.

A numerical example involving a detailed engineering computational model was presented to demonstrate the efficiency and accuracy of the proposed method. The relevant SLE analysis calculation program was compiled for the computational model above, as well as a finite element model established with ANSYS software for comparison and verification of the proposed method. The results show that the proposed method exhibited a better performance in terms of accuracy and efficiency, 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 efficiency. Moreover, the engineering example demonstrates that the SLE has great influence 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 effect 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.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the Open Fund of Key Lab of Highway Construction and Maintenance Technology in Loess Region of Shanxi Transportation Research Institute (KLTLR-Y12-9). This support is appreciated.