Improved finite beam element method for analyzing the flexural natural vibration of thin-walled box girders

Abstract

This study establishes the improved element stiffness and mass matrices of the thin-walled box girder, using a cubic Hermite polynomial shape function and based on an improved displacement function for shear-lag warping meeting the axial equilibrium condition of shear-lag warping stress and the consistency requirements of the displacement function. The improved thin-walled box girder element comprehensively considers multiple influencing factors, such as the thin-walled box girder shear-lag effect, shear deformation, and rotational inertia. The element shape function has first-order continuity at the element interface and satisfies the calculation precision with relatively few degrees of freedom. A finite beam element method program that can be used to calculate the natural vibration frequencies of thin-walled box girders was compiled based on the improved thin-walled box girder element. This program was used to calculate the natural vibration frequencies of many thin-walled box girder samples with different span-width ratios, span-height ratios, and boundary conditions. After comprehensively considering the influencing factors, the calculation results obtained from finite beam element method are in good agreement with those obtained from an ANSYS finite element calculation, thus demonstrating the rationality and validity of the improved thin-walled box girder element stiffness matrix and mass matrix proposed in this study. Finally, factors influencing the thin-walled box girder shear-lag effect and shear deformation were analyzed, providing relevant reference for project designers.

Introduction

Because shear strain exists in the top, cantilever, and bottom slabs of thin-walled box girders (TWBGs), the longitudinal strain of the portion of these slabs far from the web is less than that of the portion near the web. This phenomenon is known as the shear-lag effect and results in a nonuniform strain distribution across the transverse direction of the top, cantilever, and bottom slabs of a TWBG.^1–8 Additionally, neglecting shear deformation, as in Euler–Bernoulli beams (i.e. beams with infinite shear stiffness), might be inappropriate in some cases, for example, for TWBGs with a reduced span-to-depth ratio.^9–16 Therefore, the natural vibration characteristics of TWBGs suffer from the coupled effects of shear lag and shear deformation.

Shear-lag and shear-deformation effects in TWBGs have been popular research topics for decades. Reissner ¹⁷ first analyzed the shear-lag effect in TWBGs based on the principle of minimum potential energy. The TWBGs in his paper had a single-cell cross-section without side cantilever slabs, and a quadratic parabola was adopted as the displacement function for shear-lag warping. Reissner’s method was later extended to analyze the shear lag of TWBGs with cantilever slabs by Kuzmanovic and Graham, ¹⁸ Chang and Zheng, ¹⁹ and Chang and Yun, ²⁰ among whom Kuzmanovic and Graham ¹⁸ and Chang and Yun ²⁰ adopted quadratic and quartic parabolas, respectively. However, in most of these papers, the warping displacement functions for shear lag of the TWBG with cantilever slabs have some deficiencies because the corresponding warping stresses cannot satisfy the axial equilibrium condition. Moreover, in existing analytical methods, the shear-lag deformation state is analyzed together with the flexural deformation state of the corresponding elementary beam, which complicates shear-lag analysis. Additionally, it is inconvenient to distinguish between the shear-lag warping stress and the flexural stress of the corresponding elementary beam. Zhang and colleagues^2,3,7 established an improved displacement function for shear-lag warping in a box girder with cantilever slabs through the axial equilibrium condition for shear-lag warping stress. However, these models failed to consider shear deformation. Xu and Wu developed a new plane-stress model using the state space method for composite beams with flexible shear connections. The results show that the one-dimensional (1D) theory neglecting shear deformation underestimates the deflection of composite beams. ²¹ Berczyński and Wróblewski^9,10 formulated three analytical models describing the dynamic behavior of the composite beams: two based on Euler beam theory and one on Timoshenko beam theory. It was found out that the results obtained based on the Timoshenko beam theory model achieved the highest conformity with experimental results, both for higher and lower modes of flexural vibrations of the beam. Dilena and Morassi^22–24 further presented a refined model of composite beams considering the shear deformation and compared the results of experimental tests with those of Euler–Bernoulli and Timoshenko models. They found that the percentage errors of the Timoshenko model were less than half of those corresponding to the Euler–Bernoulli model.

In the numerical analysis of natural vibration in TWBGs, many scholars preferred to apply the finite beam element method (FBEM).^25–33 However, few scholars have comprehensively considered shear lag and shear deformation of TWBGs. In this article, the additional deflection induced by the shear-lag effect is adopted as the generalized displacement to describe the shear-lag deformation state of a TWBG and the shear-lag deformation state is analyzed as an independent fundamental deformation state. Based on the concept of generalized force corresponding to the generalized displacement for shear lag and the consistency requirements of displacement functions, improved TWBG element stiffness and mass matrices with a cubic Hermite polynomial shape function were derived. The improved TWBG element comprehensively considered multiple influencing factors, such as the shear-lag effect, shear deformation, and rotational inertia, and satisfied the calculation precision with a relatively small number of degrees of freedom, providing a foundation to improve the computational efficiency of thin-walled beam structures. Second, an FBEM program that can be used to calculate the natural vibration frequencies of TWBGs commonly used in engineering was compiled based on the improved TWBG element proposed in the study. Finally, the natural vibration frequencies of many TWBG samples with different span-width ratios, span-height ratios, and boundary conditions were calculated, and factors influencing the TWBG shear-lag effect and shear deformation were analyzed.

Displacement and strain model of TWBG

Displacement and strain models of TWBGs

Set $g_i(x,y)(i=1,2,3)$ and $g_4(x,z)$ to the displacement function for shear-lag warping of the top slab, cantilever slab, bottom slab, and web, respectively. Since the transverse shear strains at the edge of cantilever slab, the centroid of top slab, and the centroid of bottom slab are all zero, the following equations can be obtained

$$\frac{\partial g_1(x,0)}{\partial y}=\frac{\partial g_2(x,b_1+b_2)}{\partial y}=\frac{\partial g_3(x,0)}{\partial y}=0$$

(1)

where $2b_1$ , $b_2$ , and $2b_3$ denote the widths of top, cantilever, and bottom slabs, respectively, as shown in Figure 1.

Figure 1.

Dimensions of TWBG sections.

According to the displacement compatibility at the intersection of the web and top flange, the following equations are obtained

$$g_1(x,b_1)=g_2(x,b_1)=g_4(x,z_0),g_3(x,b_3)=g_4(x,z_3)$$

(2)

where $z_0$ denotes the z-coordinate of the intersections of the concrete slab and web and $z_3$ denotes the z-coordinate of the bottom slab.

Although most scholars preferred the cubic parabola,^34–37 Zhang and colleagues^2,3,7 demonstrated that a quadratic parabola was reasonable as the displacement function for shear-lag warping by investigating the in-plane shear deformation of flanges. Based on the previously mentioned boundary conditions, the displacement function for shear-lag warping of the top slab, cantilever slab, bottom slab, and web can be approximated by four quadratic parabolic equations as

$$g_i(x,y)={\psi }_i\left(y\right)U\left(x\right)i=1,2,3,4$$

(3)

$${\psi }_i\left(y\right)={\alpha }_i(y^2/b_i^2-1)+di=1,2,3,4$$

(4)

where the y of ${\psi }_2(y)$ should be replaced by $\overline{y}=b_1+b_2-\left|y\right|$ and $U(x)$ denotes the amplitude function of the displacement function for shear-lag warping. As indicated in equation (4), the constant d is added to the entire cross-section, through which the displacement function for shear-lag warping can satisfy the axial equilibrium condition.

The longitudinal displacement of the web is composited by cross-section rotation and warping. According to the beam theories in open literature, the nonlinear warping of the box girder is concentrated on top slab, bottom slab, and the cantilever slab, while the nonlinear warping of the web is insignificant. Thereby, the warping of the web can be assumed to be constant. The ${\alpha }_i(i=1,2,3,4)$ may be approximately taken as ¹²

$${\alpha }_1=1,{\alpha }_2=b_2^2/b_1^2,{\alpha }_3=z_3{b}_3^2/\left(b_1^2{z}_1\right),{\alpha }_4=0$$

(5)

where $z_1$ denotes the z-coordinates of the centroids of the top slab.

Per the axial equilibrium condition, the displacement function for shear-lag warping should meet the following equation

$$\sum _{i=1}^4\underset{A_{\mathrm{i}}}{\int }{\psi }_{\mathrm{i}}U\prime \mathrm{d}A=0$$

(6)

Substituting equation (4) into equation (6) leads to

$$d=2({\alpha }_1{A}_1+{\alpha }_2{A}_2+{\alpha }_3{A}_3)/\left(3A_0\right)$$

(7)

where $A_1=2b_1{t}_1$ , $A_2=2b_2{t}_2$ , $A_3=2b_3{t}_3$ , $A_w=A_4=2b_4{t}_4$ , $A_0=A_1+A_2+A_3+A_4$ ; $t_i(i=1,2,3,4)$ denotes the thicknesses of the top slab, cantilever slab, bottom slab, and web, respectively.

The longitudinal displacement at any point of the cross-section of the TWBG can be expressed as ³⁸

$$u_i(x,y,z)=-z\theta \left(x\right)+g_i(x,y)i=1,2,3,4$$

(8)

where $\theta (x)=w\prime (x)-{\gamma }_{\mathrm{w}}(x)$ denotes the cross-section rotation angle and $w(x)$ denotes the vertical deflection of the TWBG.

Strain model of TWBG

Based on these displacement models, the strain at any point of the TWBG’s cross-section follows

$${\epsilon }_{\mathrm{xi}}=\frac{\partial u_i}{\partial x}=-z\theta \prime +{\psi }_{i}U\prime i=1,2,3,4$$

(9)

$${\gamma }_{\mathrm{xyi}}=\frac{\partial u_i}{\partial y}=\psi {\prime }_{i}U(i=1,2,3),{\gamma }_{\mathrm{w}}=w\prime -\theta$$

(10)

where ${\epsilon }_{\mathrm{xi}}(i=1,2,3,4)$ denotes the longitudinal strain at any point of the top slab, cantilever slab, bottom slab, and web, respectively; ${\gamma }_{\mathrm{xyi}}(i=1,2,3)$ denotes the shear strain at any point of the top, cantilever, and bottom slabs, respectively; and ${\gamma }_{\mathrm{w}}$ denotes the shear strain at any point of the web.

Given these strain models, the stresses at any point of the cross-section of the TWBG can be given by

$${\sigma }_{\mathrm{xi}}=E(-z\theta \prime +{\psi }_{i}U\prime )i=1,2,3,4$$

(11)

$${\tau }_{\mathrm{xyi}}=G\psi {\prime }_{i}U(i=1,2,3){\tau }_{\mathrm{w}}=G{\gamma }_{\mathrm{w}}$$

(12)

where ${\sigma }_{\mathrm{xi}}(i=1,2,3,4)$ denotes the longitudinal stress at any point of the top slab, cantilever slab, bottom slab, and web, respectively; ${\tau }_{\mathrm{xyi}}(i=1,2,3)$ denotes the shear stress at any point on the top, cantilever, and bottom slabs, respectively; and ${\tau }_{\mathrm{w}}$ denotes the shear stress at any point of the web.

FBEM of TWBGs

TWBG element stiffness and mass matrices considering shear deformation and shear lag

According to the consistency requirements of displacement functions, the contributions of different displacement functions to the same strain function should have the same polynomial degree after finite element approximation; otherwise, false geometrical constraints will be introduced, thereby reducing the accuracy of the element and even causing shear locking. By investigating the highest derivatives of the different displacement functions in equations (9) and (10), it was found that the highest derivatives of $U(x)$ and $\theta (x)$ belong to first-order derivatives. Therefore, its shape function need only satisfy the ${\mathrm{C}}^0$ continuity condition. As indicated by equation (10), the shape function of the deflection function $w(x)$ should satisfy ${\mathrm{C}}^1$ continuity. According to the numerical calculation principle, Hermite polynomials that meet the ${\mathrm{C}}^1$ continuity condition are at least cubic. As indicated by equation (10), the contribution of the deflection function $w(x)$ to the cross-section rotation angle function $\theta (x)$ is a quadratic polynomial. Per the consistency requirements, displacement functions $U(x)$ and $\theta (x)$ should apply the quadratic polynomial to approximate. It is inappropriate to make $U(x)$ and $w(x)$ be interpolated linearly, which cannot represent the real displacement of the beam element and even cause shear locking. Therefore, elements that meet the consistency requirements need at least 10 nodal degrees of freedom, as shown in Figure 2. The element nodal displacement vector can be given by

$$\begin{array}{l}{q}_{\mathrm{u}}={\left[U_1,U_2,U_3\right]}^{\mathrm{T}},q_{\theta }={\left[{\theta }_1,{\theta }_2,{\theta }_3\right]}^{\mathrm{T}},\\ q_{\mathrm{w}}={\left[w_1,w_{1^{\prime }},w_2,w_{2^{\prime }}\right]}^{\mathrm{T}}\end{array}$$

(13)

Figure 2.

Nodal degrees of freedom of the TWBG element.

The Hermite polynomial shape function that realizes the element-edge ${\mathrm{C}}^1$ continuity condition can be expressed as

$${\phi }_{\mathrm{u}}={\phi }_{\theta }=[1-3k+2k^2,4k-4k^2,-k+2k^2]$$

(14)

$${\phi }_{\mathrm{w}}=\left[1-3k^2+2k^3,\left(k-2k^2+k^3\right)l,3k^2-2k^3,l\left(-k^2+k^3\right)\right]$$

(15)

where $k=x/l$ and l denotes element length.

$U(x)$ , $w(x)$ , and $\theta (x)$ may be approximated by the Hermite polynomial inside the element as

$$U={\phi }_{\mathrm{u}}{\mathbf{q}}_{\mathrm{u}},\theta ={\phi }_{\theta }{\mathbf{q}}_{\theta },w={\phi }_{\mathrm{w}}{\mathbf{q}}_{\mathrm{w}}$$

(16)

Substituting equation (16) into equations (9) and (10) leads to

$${\epsilon }_{\mathrm{xi}}=-z\phi {\prime }_{\theta }{\mathbf{q}}_{\theta }+{\psi }_i\phi {\prime }_{\mathrm{u}}{\mathbf{q}}_{\mathrm{u}}i=1,2,3,4$$

(17)

$${\gamma }_{\mathrm{xyi}}=\psi {\prime }_i{\phi }_{\mathrm{u}}{\mathbf{q}}_{\mathrm{u}},{\gamma }_{\mathrm{w}}=\phi {\prime }_{\mathrm{w}}{\mathbf{q}}_{\mathrm{w}}-{\phi }_{\theta }{\mathbf{q}}_{\theta }$$

(18)

Substituting equation (16) into equations (11) and (12) leads to

$${\sigma }_{\mathrm{xi}}=E(-z{\phi \prime }_{\theta }{\mathbf{q}}_{\theta }+{\psi }_i{\phi \prime }_{\mathrm{u}}{\mathbf{q}}_{\mathrm{u}})i=1,2,3,4$$

(19)

$${\tau }_{\mathrm{xyi}}=G\psi {\prime }_i{\phi }_{\mathrm{u}}{\mathbf{q}}_{\mathrm{u}},{\tau }_{\mathrm{w}}=G[{\phi \prime }_{\mathrm{w}}{\mathbf{q}}_{\mathrm{w}}-{\phi }_{\theta }{\mathbf{q}}_{\theta }]$$

(20)

The kinetic energy of the TWBG can be given by

$$T=\frac{1}{2}\underset{L}{\int }(\sum _{i=1}^4\underset{A_i}{\int }\rho {\stackrel{·}{u}}_i^2\mathrm{d}A+m{\stackrel{·}{w}}^2)\mathrm{d}x$$

(21)

where L denotes effective length of the TWBG.

The strain energy of the TWBG can be expressed as

$$V=\frac{1}{2}\underset{L}{\int }[\sum _{i=1}^4\underset{A_i}{\int }({\sigma }_{\mathrm{xi}}{\epsilon }_{\mathrm{xi}}+{\tau }_{\mathrm{xyi}}{\gamma }_{\mathrm{xyi}})\mathrm{d}A+\underset{A_{\mathrm{w}}}{\int }{\tau }_{\mathrm{w}}{\gamma }_{\mathrm{w}}\mathrm{d}A]\mathrm{d}x$$

(22)

where $m=\rho A_0$ and $\rho$ denotes the density of the TWBG.

According to equations (21) and (22), the first-order variations of the strain energy and kinetic energy can, respectively, be expressed as

$$\underset{t_1}{\overset{t_2}{\int }}\delta T\mathrm{d}t=-\underset{t_1}{\overset{t_2}{\int }}\underset{L}{\int }(m\stackrel{··}{w}\delta w+\sum _{i=1}^4\underset{A_i}{\int }{\rho }_{\mathrm{i}}{\stackrel{··}{u}}_{\mathrm{i}}\delta u_{\mathrm{i}}\mathrm{d}A)\mathrm{d}x\mathrm{d}t$$

(23)

$$\begin{array}{l}{\displaystyle \underset{t_1}{\overset{t_2}{\int }}\delta V\mathrm{d}t}=\\ {\displaystyle \underset{t_1}{\overset{t_2}{\int }}{\displaystyle \underset{L}{\int }\left[{\displaystyle \sum _{i=1}^4{\displaystyle \underset{A_i}{\int }\left({\sigma }_{xi}\delta {\epsilon }_{xi}+{\tau }_{xyi}\delta {\gamma }_{xyi}\right)\mathrm{d}A}}+{\displaystyle \underset{A_{\mathrm{w}}}{\int }{\tau }_{\mathrm{w}}\delta {\gamma }_{\mathrm{w}}\mathrm{d}A}\right]\mathrm{d}x}}\mathrm{d}t\end{array}$$

(24)

Substituting equation (16) into equation (23) leads to the TWBG element mass matrix ²⁵

$${\mathbf{M}}_{\mathrm{e}}=\left[\begin{array}{ccc}{\mathbf{M}}_{\mathrm{uu}}\hfill & {\mathbf{M}}_{\mathrm{u}\mathrm{\theta }}\hfill & {\mathbf{M}}_{\mathrm{uw}}\hfill \\ \multicolumn{1}{c}{{\mathbf{M}}_{\mathrm{\theta }\mathrm{u}}}& {\mathbf{M}}_{\mathrm{\theta }\mathrm{\theta }}\hfill & {\mathbf{M}}_{\mathrm{\theta }\mathrm{w}}\hfill \\ \multicolumn{1}{c}{{\mathbf{M}}_{\mathrm{wu}}}& {\mathbf{M}}_{\mathrm{w}\mathrm{\theta }}\hfill & {\mathbf{M}}_{\mathrm{ww}}\hfill \end{array}\right]$$

(25)

$$\begin{array}{l}{M}_{\mathrm{uu}}={\displaystyle \underset{l}{\int }{\displaystyle \sum _{i=1}^4{\displaystyle \underset{A_i}{\int }\rho {\psi }_i^2{\phi }_{\mathrm{u}}^{\mathrm{T}}{\phi }_{\mathrm{u}}\mathrm{d}A}}\mathrm{d}x},\\ M_{\mathrm{u}\theta }=-{\displaystyle \underset{l}{\int }{\displaystyle \sum _{i=1}^4{\displaystyle \underset{A_i}{\int }\rho z{\psi }_i{\phi }_{\mathrm{u}}^{\mathrm{T}}{\phi }_{\theta }\mathrm{d}A}}}\mathrm{d}x\end{array}$$

(26)

$$\begin{array}{l}{M}_{\theta \theta }={\displaystyle \underset{l}{\int }{\displaystyle \sum _{i=1}^4{\displaystyle \underset{A_i}{\int }\rho z^2{\phi }_{\theta }^{\mathrm{T}}{\phi }_{\theta }\mathrm{d}A}}\mathrm{d}x},\\ M_{\mathrm{ww}}={\displaystyle \underset{l}{\int }{\displaystyle \sum _{i=1}^4{\displaystyle \underset{A_i}{\int }\rho {\phi }_{\mathrm{w}}^{\mathrm{T}}{\phi }_{\mathrm{w}}\mathrm{d}A}}\mathrm{d}x}\end{array}$$

(27)

$$\begin{array}{l}{M}_{\mathrm{uw}}=0,M_{\theta \mathrm{w}}=0,M_{\theta \mathrm{u}}=M_{\mathrm{u}\theta }^{\mathrm{T}},\\ M_{\mathrm{w}\theta }=M_{\theta \mathrm{w}}^{\mathrm{T}},M_{\mathrm{wu}}=M_{\mathrm{uw}}^{\mathrm{T}}\end{array}$$

(28)

Substituting equations (17)–(20) into equation (24) leads to the TWBG element stiffness matrix

$${\mathbf{K}}_{\mathrm{e}}=\left[\begin{array}{ccc}{\mathbf{K}}_{\mathrm{uu}}\hfill & {\mathbf{K}}_{\mathrm{u}\mathrm{\theta }}\hfill & {\mathbf{K}}_{\mathrm{uw}}\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{\mathrm{\theta }\mathrm{u}}}& {\mathbf{K}}_{\mathrm{\theta }\mathrm{\theta }}\hfill & {\mathbf{K}}_{\mathrm{\theta }\mathrm{w}}\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_{\mathrm{wu}}}& {\mathbf{K}}_{\mathrm{w}\mathrm{\theta }}\hfill & {\mathbf{K}}_{\mathrm{ww}}\hfill \end{array}\right]$$

(29)

$${\mathbf{K}}_{\mathrm{\theta }\mathrm{w}}=-\underset{L}{\int }\underset{A_{\mathrm{w}}}{\int }G{\phi }_{\mathrm{\theta }}^{\mathrm{T}}{\phi }_{\mathrm{w}}\prime \mathrm{d}A\mathrm{d}x,{\mathbf{K}}_{\mathrm{ww}}=\underset{L}{\int }\underset{A_{\mathrm{w}}}{\int }G{{\phi }_{\mathrm{w}}\prime }^T{\phi }_{\mathrm{w}}\prime \mathrm{d}A\mathrm{d}x$$

(30)

$$\begin{array}{l}{K}_{\mathrm{uu}}={\displaystyle \underset{L}{\int }{\displaystyle \sum _{i=1}^4{\displaystyle \underset{A_i}{\int }\left(E{\psi }_i^2{{\phi }_{{\mathrm{u}}^{\prime }}}^{\mathrm{T}}{\phi }_{u^{\prime }}+G{{{\psi }^{\prime }}_i}^2{\phi }_{\mathrm{u}}^{\mathrm{T}}{\phi }_u\right)}\mathrm{d}A}}\mathrm{d}x,\\ K_{\mathrm{w}\theta }=K_{\mathrm{w}\theta }^{\mathrm{T}}\end{array}$$

(31)

$${\mathbf{K}}_{\mathrm{u}\mathrm{\theta }}=-\underset{L}{\int }\sum _{i=1}^4\underset{A_i}{\int }\mathrm{Ez}{\psi }_i{{\phi }_u\prime }^{\mathrm{T}}{\phi }_{\theta }\prime \mathrm{d}A\mathrm{d}x,{\mathbf{K}}_{\mathrm{uw}}=0,{\mathbf{K}}_{\mathrm{\theta }\mathrm{u}}={\mathbf{K}}_{\mathrm{u}\mathrm{\theta }}^{\mathrm{T}}$$

(32)

$$\begin{array}{l}{K}_{\theta \theta }={\displaystyle \underset{L}{\int }\left({\displaystyle \sum _{i=1}^4{\displaystyle \underset{A_i}{\int }{z}^{2}E{{\phi }_{{\theta }^{\prime }}}^{\mathrm{T}}{\phi }_{{\theta }^{\prime }}\mathrm{d}A}}+{\displaystyle \underset{A_{\mathrm{w}}}{\int }G{\phi }_{\theta }^{\mathrm{T}}{\phi }_{\theta }\mathrm{d}A}\right)}\mathrm{d}x,\\ K_{\mathrm{wu}}=K_{\mathrm{uw}}^{\mathrm{T}}\end{array}$$

(33)

Calculation of the natural vibration frequency of a TWBG

Based on the TWBG element stiffness matrix ${\mathbf{K}}_{\mathrm{e}}$ , element mass matrix ${\mathbf{M}}_{\mathrm{e}}$ , and element nodal displacement vector ${\mathbf{q}}_{\mathrm{e}}=\{{\mathbf{q}}_{\mathrm{u}}^{\mathrm{T}},{\mathbf{q}}_{\mathrm{\theta }}^{\mathrm{T}},{\mathbf{q}}_{\mathrm{w}}^{\mathrm{T}}{\}}^{\mathrm{T}}$ proposed by this study, the total stiffness matrix K, total mass matrix M, and the total freedom vector $\mathbf{q}$ of the TWBG can be obtained using principle of “set in the right position.” The boundary conditions commonly used by TWBGs can be expressed as follows:^15,16

1. Boundary conditions of the clamped beam

   $${U|}_{x=0,L}={w|}_{x=0,L}={\theta |}_{x=0,L}=0$$

   (34)

2. Boundary conditions of the simply supported beam

   $${U\prime |}_{x=0,L}={w|}_{x=0,L}={\theta \prime |}_{x=0,L}=0$$

   (35)

The $U\prime$ and $\theta \prime$ correspond to generalized bending moments caused by the warping deformation and bending deformation of the box girder, respectively. Therefore, the $U\prime$ and $\theta \prime$ are zero when the ends are simply supported.

Supposing the transformational relation between the TWBG total freedom vector before and after the boundary constraint as

$$\mathbf{q}=\mathbf{S}\overline{\mathbf{q}}$$

(36)

where S is the transfer matrix of the total freedom vector, which could be obtained per the boundary conditions, and $\overline{\mathbf{q}}$ is the TWBG total freedom vector after the boundary constraint.

The free vibration equation with multiple degrees of freedom can be expressed as

$$\mathbf{M}\stackrel{··}{\mathbf{u}}+\mathbf{Ku}=0$$

(37)

Assuming that

$$\mathbf{u}=\sin (\omega t+\varphi )\mathbf{q}$$

(38)

where $\omega$ denotes the structure’s natural vibration frequency and $\varphi$ denotes the initial phase.

Substituting equation (38) into equation (37) leads to

$$(-{\omega }^2\mathbf{M}+\mathbf{K})\mathbf{q}=0$$

(39)

Substituting equation (36) into equation (39) leads to

$$(-{\omega }^2\overline{\mathbf{M}}+\overline{\mathbf{K}})\overline{\mathbf{q}}=0$$

(40)

where

$$\overline{\mathbf{M}}={\mathbf{S}}^{\mathrm{T}}\mathbf{MS}\overline{\mathbf{K}}={\mathbf{S}}^{\mathrm{T}}\mathbf{KS}$$

(41)

As indicated by equation (40), the TWBG frequency equation can be given as

$$|-{\omega }^2\overline{\mathbf{M}}+\overline{\mathbf{K}}|=0$$

(42)

MATLAB was used to compile the above FBEM program for TWBGs; thus, the TWBG natural vibration frequency could be comprehensively obtained accounting for shear lag and shear deformation.

Degradation of the TWBG element stiffness matrix

• 1. When the shear-lag effect was not taken into consideration

$${\mathbf{q}}_{\mathrm{u}}={[U_1,U_2,U_3]}^{\mathrm{T}}=0$$

(43)

After the rows and columns related to the warping displacement vector ${\mathbf{q}}_{\mathrm{u}}$ were removed from the original element stiffness matrix ${\mathbf{K}}_{\mathrm{e}}$ and mass matrix ${\mathbf{M}}_{\mathrm{e}}$ , the new stiffness and mass matrices (i.e. the element stiffness matrix ${\tilde{\mathbf{K}}}_{\mathrm{e}}$ and the element mass matrix ${\tilde{\mathbf{M}}}_{\mathrm{e}}$ ), which do not consider the shear-lag effect, were formed

$${\tilde{\mathbf{K}}}_{\mathrm{e}}={\mathbf{S}}_{\mathrm{l}}{\mathbf{K}}_{\mathrm{e}}{\mathbf{S}}_{\mathrm{l}}^{\mathrm{T}},{\tilde{\mathbf{M}}}_{\mathrm{e}}={\mathbf{S}}_{\mathrm{l}}{\mathbf{M}}_{\mathrm{e}}{\mathbf{S}}_{\mathrm{l}}^{\mathrm{T}}$$

(44)

$${\left[S_{\mathrm{l}}\right]}_{7\times 10}\left(i,j+3\right)={\delta }_{ij}\left(i=1,2,\dots ,7;j=-2,-1,\dots ,7\right)$$

(45)

• 2. When shear deformation was not taken into consideration

$${\gamma }_{\mathrm{w}}=w\prime -\theta =0\Rightarrow w\prime =\theta \Rightarrow \phi {\prime }_{\mathrm{w}}{\mathbf{q}}_{\mathrm{w}}={\phi }_{\mathrm{\theta }}{\mathbf{q}}_{\mathrm{\theta }}$$

(46)

$$w\prime {|}_{x=0}=\phi {\prime }_w{\mathbf{q}}_{\mathrm{w}}{|}_{x=0}=w{\prime }_1,w\prime {|}_{x=l}=\phi {\prime }_w{\mathbf{q}}_{\mathrm{w}}{|}_{x=l}=w{\prime }_2$$

(47)

$$w\prime {|}_{x=0.5l}=\phi {\prime }_{\mathrm{w}}{\mathbf{q}}_{\mathrm{w}}{|}_{x=0.5l}=[-\frac{3}{2l},-\frac{1}{4},\frac{3}{2l},-\frac{1}{4}]{\mathbf{q}}_{\mathrm{w}}$$

(48)

$$\begin{array}{l}\theta {|}_{x=0}={\phi }_{\theta }{q}_{\theta }{|}_{x=0}={\theta }_1,\\ \theta {|}_{x=0.5l}={\phi }_{\theta }{q}_{\theta }{|}_{x=0.5l}={\theta }_2,{\theta }_{x=l}={\phi }_{\theta }{q}_{\theta }{|}_{x=l}={\theta }_3\end{array}$$

(49)

As indicated by equations (46)–(49), the following could be obtained

$${\theta }_1=w{\prime }_1,{\theta }_2=[-\frac{3}{2l},-\frac{1}{4},\frac{3}{2l},-\frac{1}{4}]{\mathbf{q}}_{\mathrm{w}},{\theta }_3=w{\prime }_2$$

(50)

As indicated by equation (50), the element stiffness matrix ${\tilde{\mathbf{K}}}_{\mathrm{e}}$ and the element mass matrix ${\tilde{\mathbf{M}}}_{\mathrm{e}}$ that do not consider shear deformation were

$${\stackrel{⌢}{K}}_{\mathrm{e}}=S_{\mathrm{d}}{K}_{\mathrm{e}}{S}_{\mathrm{d}}^{\mathrm{T}},{\stackrel{⌢}{M}}_e=S_{\mathrm{d}}{M}_{\mathrm{e}}{S}_{\mathrm{d}}^{\mathrm{T}}$$

(51)

where ${\mathbf{S}}_{\mathrm{d}}$ is the $7\times 10$ transfer matrix that was obtained by equation (50).

After obtaining the TWBG element stiffness and mass matrices not considering the shear-lag effect and shear deformation, the method described in section “Calculation of the natural vibration frequency of a TWBG” was used to calculate the TWBG natural vibration frequency without considering shear lag or shear deformation.

Analysis of calculation examples

The finite element method (FEM) was carried out by utilizing the finite element program, ANSYS, which was used to verify the FBEM’s results. The top slab, bottom slab, and web of the box beam were simulated by employing SHELL43 shell elements. The boundary conditions of the end of the finite element models were simulated by constraining the degrees of freedom in both the vertical and transverse directions if the ends were simply supported, and in the vertical, transverse, and longitudinal directions if the ends were clamp-supported.

In the first place, a numerical example by applying Carrera Unified Formulation (CUF) ³⁹ was quoted to compare with the FBEM’s result. This study focuses on a particular class of CUF models that makes use of Lagrange polynomials to discretize cross-sectional displacement variables. This class of models is referred to as component-wise (CW) in recent works. According to the CW approach, each structural component (e.g. columns, walls, frame members, and floors) can be modeled by means of the same 1D formulation. In this study, the hollow-square section beam was investigated. The dimension of the sides were $b_1=b_3=0.1\mathrm{m}$ , $b_2=0$ , $h=0.2\mathrm{m}$ , $t_2=\beta =0$ , and $t_1=t_3=0.02\mathrm{m}$ as shown in Figure 1. The beam was made of an aluminum alloy with the following characteristics: $E=75\mathrm{GPa}$ , $\upsilon =0.3$ , and $\rho =2.7\times {10}^3\mathrm{kg}/{\mathrm{m}}^3$ . The beam’s boundary condition is clamped-free. The first four flexural natural frequencies were analyzed in Table 1. The result of FBEM model was shown in the second columns. Two different CW models applying CUF were given in the third and fourth columns, while FEM models performed both by MSC Nastran and ANSYS were shown in the fifth and sixth columns. For each model, the number of degree of freedoms (DOFs) was given in the third row of Table 1.

Table 1.

Comparison of calculation results between the FBEM and CW.

DOFs	FBEM	CW models		FEM models
	FBEM Model (121)	Model A (4680)	Model B (6840)	NAS2D (40,000)	ANSYS (12,280)
Model 1 f	53.2	53.7	53.6	53.5	53.2
Model 2 f	293.7	302.2	301.3	298.2	298.2
Model 3 f	711.8	741.9	735.9	723.5	723.7
Model 4 f	1198.9	1251.7	1229.8	1193.1	1193.9

FBEM: finite beam element method; CW: component-wise; FEM: finite element method; DOFs: degree of freedoms.

“f” denotes the flexural mode.

Table 1 shows that the FBEM’s first four flexural natural frequencies were in good agreement with the CW’s results and FEM programs’ results. Meanwhile, the FBEM’s numerical cost is relatively small to compute the box beam’s flexural natural frequencies.

Furthermore, in order to verify the accuracy of the FBEM in section “FBEM of TWBGs,” the natural frequencies of two groups of simply supported TWBGs (TWBG-1 and TWBG-2) and two groups of clamp-supported TWBGs (TWBG-3 and TWBG-4) were calculated using both the FEM and FBEM. Five span-width ratios and five span-height ratios were selected for each group of TWBGs. The mechanical and geometrical parameters of the TWBG samples were as follows: $E=3.5\times {10}^{10}\mathrm{N}/{\mathrm{m}}^2$ , $\mu =0.18$ , $\rho =2400\mathrm{kg}/{\mathrm{m}}^3$ , $b_1=b_3$ , $b_2=0.1\mathrm{m}$ , $t_1=t_2=t_3=0.025\mathrm{m}$ , $t_4=0.04\mathrm{m}$ , and $L_1=L_2=L_3=L_4=7\mathrm{m}$ , where $\mu$ denotes Poisson’s ratio of the material, $r_{\mathrm{w}}=L/2b_1$ , and $r_{\mathrm{h}}=L/h$ denote the span-width ratio and the span-height ratio of the TWBG, respectively.

Tables 2–5 contain comparisons of the calculation results for natural frequencies using the FBEM and ANSYS models. $R_{\mathrm{AN}}$ denotes the calculation results obtained from ANSYS, $R_{\mathrm{FB}}$ is the calculation results obtained from the FBEM with both shear deformation and the shear-lag effect considered, $R_{\mathrm{SL}}$ denotes the calculation results obtained from the FBEM without considering the shear-lag effect, and $R_{\mathrm{SD}}$ denotes the calculation results obtained from the FBEM without considering shear deformation. $e_{\mathrm{FB}}=100(R_{\mathrm{FB}}-R_{\mathrm{AN}})/R_{\mathrm{AN}}$ represents the FBEM calculation error and $e_{\mathrm{SL}}=100(R_{\mathrm{SL}}-R_{\mathrm{FB}})/R_{\mathrm{FB}}$ and $e_{\mathrm{SD}}=100(R_{\mathrm{SD}}-R_{\mathrm{FB}})/R_{\mathrm{FB}}$ denote the shear-lag and shear-deformation effects, respectively.

Table 2.

Comparison of calculation results between the FBEM and ANSYS models (TWBG-1).

$r_{\mathrm{w}}$	Computation methods	Natural frequencies (Hz)
		1st	2nd	3rd	4th	5th
4.00	$R_{\mathrm{AN}}$	31.326	73.679	126.040	183.710	243.850
	$R_{\mathrm{FB}}$	31.134	72.747	124.869	183.958	248.650
	$R_{\mathrm{SL}}$	35.819	88.319	154.545	228.556	307.259
	$R_{\mathrm{SD}}$	33.801	83.222	149.382	230.093	325.201
6.00	$R_{\mathrm{AN}}$	33.045	80.548	140.490	207.640	278.850
	$R_{\mathrm{FB}}$	33.287	80.590	140.342	208.108	281.814
	$R_{\mathrm{SL}}$	35.842	90.590	161.761	243.099	330.930
	$R_{\mathrm{SD}}$	35.669	91.115	166.281	257.550	363.448
8.00	$R_{\mathrm{AN}}$	33.584	83.915	148.790	222.520	301.560
	$R_{\mathrm{FB}}$	33.963	84.395	149.265	223.440	304.213
	$R_{\mathrm{SL}}$	35.498	91.002	164.570	249.990	343.295
	$R_{\mathrm{SD}}$	36.031	94.279	174.975	273.940	388.767
10.00	$R_{\mathrm{AN}}$	33.643	85.480	153.520	231.890	316.680
	$R_{\mathrm{FB}}$	34.063	86.160	154.350	233.130	319.314
	$R_{\mathrm{SL}}$	35.061	90.714	165.487	253.335	350.156
	$R_{\mathrm{SD}}$	35.874	95.251	178.995	282.962	404.296
12.00	$R_{\mathrm{AN}}$	33.499	86.102	156.160	237.770	326.830
	$R_{\mathrm{FB}}$	33.925	86.870	157.178	239.219	329.482
	$R_{\mathrm{SL}}$	34.613	90.141	165.494	254.825	353.992
	$R_{\mathrm{SD}}$	35.534	95.213	180.484	287.481	413.222

FBEM: finite beam element method; TWBG-1: thin-walled box girder 1.

Table 3.

Comparison of calculation results between the FBEM and ANSYS models (TWBG-2).

$r_{\mathrm{h}}$	Computation methods	Natural frequencies (Hz)
		1st	2nd	3rd	4th	5th
4.00	$R_{\mathrm{AN}}$	111.940	234.540	376.090	517.720	615.220
	$R_{\mathrm{FB}}$	112.634	237.189	382.104	529.355	627.299
	$R_{\mathrm{SL}}$	113.800	239.933	386.375	535.409	640.122
	$R_{\mathrm{SD}}$	160.162	393.076	674.874	973.286	1274.333
6.00	$R_{\mathrm{AN}}$	90.741	200.900	328.830	462.430	597.190
	$R_{\mathrm{FB}}$	90.713	201.818	332.144	470.704	614.817
	$R_{\mathrm{SL}}$	92.239	206.179	339.587	481.342	628.580
	$R_{\mathrm{SD}}$	114.540	292.769	525.294	789.735	1071.968
8.00	$R_{\mathrm{AN}}$	75.959	174.460	291.040	414.870	540.850
	$R_{\mathrm{FB}}$	75.779	174.637	292.776	420.481	554.237
	$R_{\mathrm{SL}}$	77.437	179.998	302.655	435.132	573.632
	$R_{\mathrm{SD}}$	89.980	233.534	427.039	655.322	907.850
10.00	$R_{\mathrm{AN}}$	65.217	153.650	260.320	375.250	493.260
	$R_{\mathrm{FB}}$	65.066	153.560	261.225	379.114	503.504
	$R_{\mathrm{SL}}$	66.744	159.470	272.789	396.876	527.526
	$R_{\mathrm{SD}}$	74.477	194.643	359.202	557.080	780.534
12.00	$R_{\mathrm{AN}}$	57.123	137.050	235.110	342.080	452.890
	$R_{\mathrm{FB}}$	57.055	136.936	235.606	344.817	460.838
	$R_{\mathrm{SL}}$	58.700	143.101	248.257	364.856	488.481
	$R_{\mathrm{SD}}$	63.760	167.218	310.080	483.671	682.048

FBEM: finite beam element method; TWBG-2: thin-walled box girder 2.

Table 4.

Comparison of calculation results between the FBEM and ANSYS models (TWBG-3).

$r_{\mathrm{w}}$	Computation methods	Natural frequencies (Hz)
		1st	2nd	3rd	4th	5th
4.00	$R_{\mathrm{AN}}$	16.273	56.549	108.890	167.290	228.600
	$R_{\mathrm{FB}}$	16.463	56.285	107.660	166.249	230.724
	$R_{\mathrm{SL}}$	17.247	64.286	131.177	209.019	292.104
	$R_{\mathrm{SD}}$	16.829	60.176	120.569	194.685	282.366
6.00	$R_{\mathrm{AN}}$	16.232	59.274	118.740	187.040	259.820
	$R_{\mathrm{FB}}$	16.510	59.805	119.087	187.434	261.760
	$R_{\mathrm{SL}}$	16.857	63.965	133.413	217.109	309.000
	$R_{\mathrm{SD}}$	16.783	63.226	131.864	217.250	316.870
8.00	$R_{\mathrm{AN}}$	16.019	59.988	123.280	198.080	279.120
	$R_{\mathrm{FB}}$	16.302	60.738	124.253	199.310	281.479
	$R_{\mathrm{SL}}$	16.489	63.166	133.372	219.782	316.422
	$R_{\mathrm{SD}}$	16.520	63.678	136.090	228.487	337.239
10.00	$R_{\mathrm{AN}}$	15.776	59.926	125.190	204.130	291.030
	$R_{\mathrm{FB}}$	16.046	60.730	126.418	205.767	293.685
	$R_{\mathrm{SL}}$	16.161	62.276	132.537	220.250	319.635
	$R_{\mathrm{SD}}$	16.228	63.287	137.197	233.423	348.143
12.00	$R_{\mathrm{AN}}$	15.542	59.557	125.800	207.320	298.330
	$R_{\mathrm{FB}}$	15.795	60.357	127.109	209.136	301.109
	$R_{\mathrm{SL}}$	15.871	61.408	131.418	219.713	320.731
	$R_{\mathrm{SD}}$	15.953	62.619	136.934	235.065	353.380

FBEM: finite beam element method; TWBG-3: thin-walled box girder 3.

Table 5.

Comparison of calculation results between the FBEM and ANSYS models (TWBG-4).

$r_{\mathrm{h}}$	Computation methods	Natural frequencies (Hz)
		1st	2nd	3rd	4th	5th
4.00	$R_{\mathrm{AN}}$	66.151	204.400	359.090	513.580	532.840
	$R_{\mathrm{FB}}$	65.839	205.120	363.555	525.574	544.417
	$R_{\mathrm{SL}}$	66.107	207.004	367.942	532.549	558.254
	$R_{\mathrm{SD}}$	72.333	260.677	513.900	796.007	1088.644
6.00	$R_{\mathrm{AN}}$	49.159	163.940	302.680	446.490	587.730
	$R_{\mathrm{FB}}$	48.861	163.654	304.135	452.876	603.998
	$R_{\mathrm{SL}}$	49.136	166.075	310.613	464.106	619.763
	$R_{\mathrm{SD}}$	51.766	194.026	400.041	645.962	915.394
8.00	$R_{\mathrm{AN}}$	39.292	136.580	260.510	393.090	525.830
	$R_{\mathrm{FB}}$	39.098	136.176	260.909	396.763	536.993
	$R_{\mathrm{SL}}$	39.359	138.794	268.621	411.029	557.956
	$R_{\mathrm{SD}}$	40.729	155.036	325.925	537.055	776.164
10.00	$R_{\mathrm{AN}}$	32.810	116.970	228.170	350.200	474.490
	$R_{\mathrm{FB}}$	32.718	116.686	228.257	352.503	482.524
	$R_{\mathrm{SL}}$	32.960	119.330	236.612	368.760	507.309
	$R_{\mathrm{SD}}$	33.761	129.474	274.826	457.628	668.549
12.00	$R_{\mathrm{AN}}$	28.218	102.290	202.790	315.340	431.700
	$R_{\mathrm{FB}}$	28.212	102.192	202.876	316.963	437.738
	$R_{\mathrm{SL}}$	28.437	104.783	211.500	334.433	465.192
	$R_{\mathrm{SD}}$	28.939	111.432	237.779	398.231	585.330

FBEM: finite beam element method; TWBG-4: thin-walled box girder 4.

The shear-deformation and shear-lag effects of TWBG were studied by analyzing the first five orders of the natural frequencies. Tables 2–5 and Figures 3–5 show that:

1. After considering influencing factors, such as shear deformation, shear lag, and rotational inertia, the calculation results obtained from FBEM are in good agreement with those obtained from ANSYS finite element calculation, with a calculation error $e_{\mathrm{FB}}$ less than 3%. This result demonstrates the rationality and validity of the improved TWBG element stiffness and mass matrices proposed in this study and provides a theoretical basis for applying FBEM to the dynamic analysis of TWBGs.

2. Both the shear-lag effect and shear deformation increase as the order of the natural vibration frequency of the TWBG rises. The shear-lag and shear-deformation effects of clamp-supported TWBGs at low-order frequencies are relatively small, indicating that the low-order natural vibration in clamp-supported TWBGs is mainly flexural deflection caused by the cross-section’s rotation.

3. The shear-lag effect significantly increases as the span-width ratio of the TWBG decreases. When the span-width ratio of the TWBG decreases from 12 to 4, the shear-lag effect is increased fourfold. When the span-width ratio is 4.0, the maximum shear-lag effect reaches 26.6%, indicating that the shear-lag effect of TWBGs with a relatively small span-width ratio cannot be ignored.

4. The shear-lag effect increases slightly as the span-height ratio of the TWBG increases; meanwhile, the shear-lag effect is not significant. This indicates that with a relatively small span-height ratio, the shear-lag effect is not a major concern.

5. The shear-deformation effect significantly increases as the span-height ratio of the TWBG decreases. When the span-height ratio of the TWBG decreases from 12 to 4, the shear-deformation effect is increased by 2.5 times. When the span-height ratio is 4.0, the shear-deformation effect reaches 103.15%, indicating that the shear-deformation effect of TWBG with a relatively small span-height ratio cannot be ignored.

6. While the span-width ratio varies, the shear-deformation effect increases when the mode order of natural frequency increases, indicating that the shear-deformation effect cannot be ignored with a higher mode.

Figure 3.

Relationship between the calculation errors $e_{\mathrm{FB}}$ and the mode orders of natural vibration of (a) TWBG-1, (b) TWBG-2, (c) TWBG-3, and (d) TWBG-4.

Figure 4.

Relationship between the shear-lag effect $e_{\mathrm{SL}}$ and the mode orders of natural vibration of (a) TWBG-1, (b) TWBG-2, (c) TWBG-3, and (d) TWBG-4.

Figure 5.

Relationship between the shear-deformation effect $e_{\mathrm{SD}}$ and the mode orders of natural vibration of (a) TWBG-1, (b) TWBG-2, (c) TWBG-3, and (d) TWBG-4.

Conclusion

After considering many influencing factors, such as shear deformation, shear lag, and rotational inertia, an improved TWBG’s element mass and stiffness matrices were derived using Hamilton’s principle with a cubic Hermite polynomial shape function. An FBEM program was implemented employing the proposed TWBG element and can be used to calculate the natural vibration frequencies of TWBGs. The natural vibration frequencies of many TWBG samples with different span-width ratios, span-height ratios, and boundary conditions were calculated, leading to the following conclusions:

1. The FBEM calculation results are in good agreement with the ANSYS results, thus demonstrating the rationality and validity of the improved TWBG element mass and stiffness matrices proposed in this study and providing the basis for further application of the FBEM in dynamic calculation for TWBGs.

2. Adopting the additional deflection as the generalized displacement for shear lag has very clear, physically significant meaning and is especially easy for engineers to accept in practice.

3. The shear-lag effect increases with the order of the natural vibration frequency of the TWBG and increases as the span-width ratio decreases. The influence of the shear-lag effect on high-order natural vibration frequencies of TWBGs with relatively small span-width ratios cannot be ignored.

4. The shear-deformation effect increases as the order of the natural vibration frequency of the TWBG increases and as the span-height ratio decreases. The influence of the shear-deformation effect on high-order natural vibration frequencies of TWBGs with relatively small span-height ratios cannot be ignored.

5. The shear-lag effect and shear deformation of clamp-supported TWBGs at low-order frequencies are relatively small, indicating that the low-order vibration in clamp-supported TWBGs is mainly flexural deflection caused by the cross-section’s rotation.