Finite Beam Element for Curved Steel–Concrete Composite Box Beams Considering Time-Dependent Effect

Abstract

Curved steel–concrete composite box beams are widely used in urban overpasses and ramp bridges. In contrast to straight composite beams, curved composite box beams exhibit complex mechanical behavior with bending–torsion coupling, including constrained torsion, distortion, and interfacial biaxial slip. The shear-lag effect and curvature variation in the radial direction should be taken into account when the beam is sufficiently wide. Additionally, long-term deflection has been observed in curved composite box beams due to the shrinkage and creep effects of the concrete slab. In this paper, an equilibrium equation for a theoretical model of curved composite box beams is proposed according to the virtual work principle. The finite element method is adopted to obtain the element stiffness matrix and nodal load matrix. The age-adjusted effective modulus method is introduced to address the concrete creep effects. This 26-DOF finite beam element model is able to simulate the constrained torsion, distortion, interfacial biaxial slip, shear lag, and time-dependent effects of curved composite box beams and account for curvature variation in the radial direction. An elaborate finite element model of a typical curved composite box beam is established. The correctness and applicability of the proposed finite beam element model is verified by comparing the results from the proposed beam element model to those from the elaborate finite element model. The proposed beam element model is used to analyze the long-term behavior of curved composite box beams. The analysis shows that significant changes in the displacement, stress and shear-lag coefficient occur in the curved composite beams within the first year of loading, after which the variation tendency becomes gradual. Moreover, increases in the central angle and shear connection stiffness both reduce the change rates of displacement and stress with respect to time.

1. Introduction

Curved steel–concrete composite box beams have gradually become one of the main design types for urban overpasses and ramp bridges due to their light weight, high torsional rigidity, and rapid construction. However, in contrast to straight beam bridges, typical bending–torsion coupling mechanical behavior can be found in curved beam bridges [1]. Hence, when the structure is subjected to a vertical load, such as its own weight, certain bending and torsional effects are generated, leading to torsional and distortional warping. For a steel–concrete composite structure, longitudinal and transverse slips exist at the interface between the steel beam and the concrete slab of the curved composite box beam. In addition, for large-curvature curved beams with a width-to-radius ratio greater than 1/10, the change in the radius of curvature along the radial direction of the beam and the shear-lag effect caused by the wide flange still need to be considered.

For the numerical analysis of curved beams, the use of a three-dimensional finite element model can accurately simulate its mechanical behavior. However, this modeling and analysis process is complex and requires high computer performance, and the calculation efficiency with a three-dimensional model is not as good as that with a one-dimensional model. The study of one-dimensional theoretical models of curved beams first began in the 1960s. Vlasov [2] established governing differential equations for curved beams based on the three major equations of mechanics, which laid the foundation for subsequent theoretical research. Later, some Japanese scholars [3,4,5] conducted related research on one-dimensional theoretical models for curved box beams and proposed relevant design proposals for the behavior of shear lag and distortion. At present, the main idea of the research on one-dimensional beam models of curved beams is to introduce functions with relevant degrees of freedom (DOFs) on the basis of Vlasov’s curved beam theoretical model to enable corresponding analyses.

The time-dependent (shrinkage and creep) effects of concrete have a significant influence on the mechanical behavior in steel–concrete composite beams. Due to the effects of concrete shrinkage and creep, the deformation of the structure under long-term loads is often much greater than that under short-term loads. The long-term deflection of composite beams has become a key issue for composite beam design and analysis. The authors of Reference [6] developed a two-node finite beam element with 26 DOFs for a curved composite box girder considering constrained torsion, distortion, shear lag, biaxial slip at the interface, and curvature variation along the radial direction. This paper introduces a constitutive relationship for concrete shrinkage and creep on the basis of the proposed beam finite element and further develops a new finite beam element that can consider the long-term behavior of curved composite box beams.

The shrinkage behavior of concrete is not related to its stress state; thus, shrinkage can be easily achieved by applying an initial strain to the model. The creep behavior of concrete is directly related to the stress state of concrete. Research on the long-term behavior of composite beams began in the 1970s. Bazant et al. [7] first simplified composite beams as a completely rigid interface and established a one-dimensional theoretical model of long-term behavior by using the three major equations of material mechanics. Afterward, Italian and Australian scholars became the primary researchers studying the long-term behavior of composite beams. Tarantino and Dezi used a stepwise calculation method to analyze the time-dependent effects of concrete, deduced a theoretical equation for simply supported composite beams, considering the time-dependent effects [8], and then derived theoretical equations for continuous composite beams, using the force method [9,10]. Then, Dezi et al. [11,12] adopted a single-step method (effective modulus method, average stress method, and age-adjusted effective modulus method) to consider the time-dependent effects of concrete and established a one-dimensional simplified theoretical model of the long-term behavior of simply supported and continuous composite beams. After entering the 21st century, based on the governing differential equations for composite beams that can consider the shear-lag behavior of concrete slabs, Dezi et al. [13] introduced a time-dependent constitutive model and established a one-dimensional theoretical model that can consider the shear-lag effect and the long-term behavior. Based on this model, the long-term behavior of composite beams was studied under prestressed loads [14]. With the increasing application of finite element technology in structural numerical analysis, a beam element model of composite beams was also proposed for analyzing their long-term behavior. Gara et al. [15] introduced a shape function to develop a beam element with 13 DOFs that can consider the interfacial slip in composite beams and the shear-lag effect of concrete slabs and analyzed the long-term behavior of the structure by using a stepwise calculation method. Ranzi and Bradford [16] proposed beam elements for the long-term behavior of composite beams, considering interfacial slip and concrete time-dependent effects according to the direct stiffness method and single-step method. Later, Ranzi, and Bradford [17] and Nguyen et al. [18] proposed long-term performance beam elements, considering interfacial slip and the time-dependent effects of concrete on the basis of the direct stiffness method and the stepwise calculation method. The abovementioned studies are all on beam element models that can analyze the long-term behavior of straight composite beams. Due to the complex mechanical behavior and time-dependent effect of curved composite box beams, it is challenging to develop a finite beam element for curved composite box beams that can account for constrained torsion, distortion, shear lag, interfacial biaxial slip, time-dependent effects, and curvature variation in the radial direction.

Based on the Vlasov’s one-dimensional theoretical model of curved beams, this paper introduces a torsional warping function, a distortional angle, an interfacial longitudinal and an interfacial transverse sliding function, and a shear-lag intensity function for concrete slabs and a shear-lag intensity function for steel beams and adopts the age-adjusted effective modulus method to analyze the effects of concrete creep. Moreover, this study establishes a finite beam element for curved composite box beams that can consider constrained torsion, distortion, shear lag, interfacial biaxial slip, time-dependent effects, and curvature variation in the radial direction. Finally, this study analyzes the long-term mechanical behavior of curved composite box girders by use of the developed finite beam element model.

2. Two-Node Finite Beam Element with 26 DOFs for Curved Composite Box Beams

2.1. Basic Hypotheses of the Model

The basic hypotheses made in the derivation process in this paper are as follows: (1) The sections of the concrete slab are rectangular, and the slab’s curvature is the same as that of the undeformed structure; (2) the vertical displacements of the concrete and the steel beam are the same; (3) the curvature radius of the composite beam does not change along the length of the beam; (4) the stiffness of the shear connections of the composite beam is uniformly distributed in the longitudinal and transverse directions; (5) the shear stress produced by the distortional effect is ignored; (6) shear deformation produced by bending and distortional warping is ignored; and (7) the shear-lag effect is taken into account under vertical bending.

2.2. Geometric Dimensions and Coordinate System of the Curved Composite Box Beam

The geometric dimensions and coordinate system of the curved composite box beam are shown in Figure 1, wherein oxyz is a three-dimensional flow coordinate system that passes through the centroid line, the ox-axis is parallel to the concrete slab and points to the curvature center of the box beam, the oy-axis is straight downward, and the oz-axis coincides with the undeformed axis of the beam. Moreover, the tangential directions of every thin-walled bar in the sections of the box beam are set as the s-axis, and the normal directions of every thin-walled bar in the sections of the box beam are set as the n-axis.

2.3. Displacement Modes and Strain Components of Curved Composite Box Beams

The displacements along the ox, oy, and oz axes of the centroid of the curved composite box beam are denoted u(z), v(z), and w(z), respectively. The torsion angle around the torsion center and the distortion angle around the distortion center are denoted θ(z) and θ_d(z), respectively. The transverse and longitudinal slip of the steel beam and concrete slab are denoted a_kΩ_x(z) and a_kΩ_s(z), respectively, where k is s or c, which represent the steel beam or concrete slab, respectively; a_s = 1 and a_c = −1; and Ω_x(z) and Ω_s(z) are the transverse and longitudinal slip functions, respectively. According to Ullmanski’s second torsion theory, an independent function of torsional warping displacement, β(z), is introduced. Additionally, f_k (k = s,c) is the shear-lag warping strength functions that represent the steel beam or concrete slab. Based on the coordinate system shown in Figure 1, according to the curved shell theory, combining the introduced transverse and longitudinal slip functions and the shear-lag warping strength functions of the steel beam and concrete slab, the longitudinal displacement W_zk(z,s), tangential displacement W_sk(z,s), and normal displacement W_nk(z,s) at any point of the composite box beam can be obtained as follows:

$$\begin{array}{l}{W}_{zk}(z,s)=w-(u^{\prime }+w{k}_0+a_k{\mathsf{\Omega }}_z{k}_0)x-v^{\prime }y+a_k{\mathsf{\Omega }}_z\\ -\omega (x,y)(\beta +v^{\prime }{k}_0)-{\omega }_{\mathrm{d}}(x,y){\theta }_{\mathrm{d}}{}^{\prime }+{\psi }_k(x)f_k\end{array}$$

(1)

$$W_{sk}(z,s)=-(u+a_k{\mathsf{\Omega }}_x)\sin \alpha -v\cos \alpha +\theta {\rho }_s+{\theta }_{\mathrm{d}}{D}_s$$

(2)

$$W_{nk}(z,s)=(u+a_k{\mathsf{\Omega }}_x)\cos \alpha -v\sin \alpha -\theta {\rho }_n-{\theta }_{\mathrm{d}}{D}_n$$

(3)

where k₀ = 1/R₀ is the initial curvature of the centroid line of the box beam sections; R₀ is the curvature radius of the centroid line of the box beam sections; α is the angle from the x-axis to the n-axis, which is anticlockwise positive; ω(x,y) and ω_d(x,y) are the torsional warping and the distortional warping principal polar coordinates, respectively; and ψ_k(x) is the shape function for the shear-lag warping of the concrete slab or steel beam, the expressions of which are as shown in Appendix A.

Transformed sections should be adopted for calculating geometric characteristics. For the calculation of the axial and bending characteristics of the composite beam, the thickness of the concrete slab is unchanged, whereas the width decreases. For the calculation of torsional and distortional effects of the composite beam, the width of the concrete slab remains unchanged, whereas the thickness decreases.

In addition, according to the literature [19,20], in Equations (1)–(3), ${\rho }_s=\left(y-y_{\mathrm{s}}\right)\sin \alpha -\left(x-x_{\mathrm{s}}\right)\cos \alpha$ is the tangential displacement distribution of the cross-sections of the composite beam under a unit torsion angle; ${\rho }_n=\left(y-y_{\mathrm{s}}\right)\cos \alpha +\left(x-x_{\mathrm{s}}\right)\sin \alpha$ is the normal displacement distribution of the cross-sections of the composite beam under a unit torsion angle; $D_s=\left[\left(y-y_{\mathrm{d}}\right)\sin \alpha -\left(x-x_{\mathrm{d}}\right)\cos \alpha \right]{\mathsf{\Psi }}_{\mathrm{d}}$ is the tangential displacement distribution of the composite beam under a unit distortion angle; $D_n=\left[\left(y-y_{\mathrm{d}}\right)\cos \alpha +\left(x-x_{\mathrm{d}}\right)\sin \alpha \right]{\mathsf{\Psi }}_{\mathrm{d}}$ is the normal displacement distribution of the composite beam under a unit distortion angle; (x_s, y_s) and (x_d, y_d) are the coordinates of the torsion center, S, and the distortion center, D, of the transformed sections of the composite beam, respectively; ${\mathsf{\Psi }}_{\mathrm{d}}=-1/\left(1+\nu \right)$ is the distortion parameter of the top and bottom plates, and ${\mathsf{\Psi }}_{\mathrm{d}}=\nu /\left(1+\nu \right)$ is the distortion parameter of the web, in which $\nu =\left(2h_1{b}_{\mathrm{t}}+2h_3{b}_{\mathrm{s}}\right)/\left(h_{2}c+h_{4}c\right)$ is a constant related only to the section shape [20]; h₁ and h₃ are the vertical distances from the distortion center to the top and bottom plates, respectively; and h₂ and h₄ are the vertical distances from the distortion center to the left and right webs, respectively.

Note that a_s is 1 for the steel beam and −1 for the concrete slab. According to Equations (1) and (2), the relative slips on the interface between the steel beam and concrete slab are expressed as follows:

$$u_{\mathrm{sp}}=2{\mathsf{\Omega }}_x$$

(4)

$$w_{\mathrm{sp}}=2{\mathsf{\Omega }}_z$$

(5)

where u_sp and w_sp denote the transverse and longitudinal interfacial slip, respectively.

Based on the abovementioned displacement mode, according to the curved shell theory and taking the curvature variation in the radial direction of the beam into consideration, the normal strain component at any point, P, in the three-dimensional flow coordinate system is expressed as follows:

$$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{p}k}& =\frac{\partial W_{zk}}{\rho (x)\partial \varphi }+\frac{W_s\sin \alpha }{\rho (x)}-\frac{W_n\cos \alpha }{\rho (x)}\hfill \\ & =\frac{R_0}{\rho (x)}\left(\frac{\partial W_{zk}}{R_0\partial \varphi }+\frac{W_s\sin \alpha }{R_0}-\frac{W_n\cos \alpha }{R_0}\right)\hfill \\ & =\frac{R_0}{\rho (x)}[w^{\prime }-x{u}^{″}-y{v}^{″}+a_k(1-k_{0}x){\mathsf{\Omega }}_z{}^{\prime }-k_{0}u-k_{0}x{w}^{\prime }-a_k{k}_0{\mathsf{\Omega }}_x+\hfill \\ & k_0(y-y_s)\theta -\omega ({\beta }^{\prime }+k_0{v}^{″})+k_0(y-y_d){\mathsf{\Psi }}_d{\theta }_{\mathrm{d}}-{\omega }_{\mathrm{d}}{\theta }_{\mathrm{d}}{}^{″}+{\psi }_k{f}_k{}^{\prime }]\hfill \end{array}$$

(6)

According to Hypothesis (5), based on the curved shell theory, the shear strain component at any point, P, of the cross-section is expressed as follows:

$${\gamma }_{\mathrm{p}k}=\frac{R_0}{\rho (x)}\left[({\theta }^{\prime }+v^{\prime }{k}_0)r^{*}-(\beta +v^{\prime }{k}_0)\frac{\partial \omega }{\partial s}+{\psi }_{k,x}{f}_k\right]$$

(7)

where r* denotes the vertical distance between the torsion center and the tangent line of point, P.

$$\frac{\partial \omega }{\partial s}=r^{*}-\frac{\mathsf{\Omega }}{\tilde{t}{\displaystyle \oint (\mathrm{d}s/\tilde{t})}}$$

(8)

where $\mathsf{\Omega }={\displaystyle \oint \tilde{t}\mathrm{d}s}$ is twice the surrounding area, and $\tilde{t}$ is the thickness of each thin-walled plate of the transformed section.

2.4. Equilibrium Equation for Curved Composite Box Beams

The equilibrium equation for a curved steel–concrete composite beam can be derived according to the virtual work principle as follows:

$$\begin{array}{cc}\hfill \mathsf{\delta }\mathsf{\Pi }=& {\displaystyle \underset{\mathsf{\Phi }}{\int }{\displaystyle \underset{A_{\mathrm{s}}}{∯}\mathsf{\delta }{\mathit{\epsilon }}_{\mathbf{s}}^{\mathrm{T}}{\mathit{\sigma }}_{\mathbf{s}}\mathrm{d}a}\rho (x)\mathrm{d}\varphi }+{\displaystyle \underset{\mathsf{\Phi }}{\int }{\displaystyle \underset{A_{\mathrm{c}}}{∯}\mathsf{\delta }{\mathit{\epsilon }}_{\mathbf{c}}^{\mathrm{T}}{\mathit{\sigma }}_{\mathbf{c}}\mathrm{d}a}\rho (x)\mathrm{d}\varphi }+{\displaystyle \underset{\mathsf{\Phi }}{\int }{\displaystyle \underset{2b_{\mathrm{ts}}}{\int }\mathsf{\delta }{\mathit{d}}_{\mathbf{slip}}^{\mathrm{T}}{\mathit{q}}_{\mathbf{sh}}\mathrm{d}x}\rho (x)\mathrm{d}\varphi }+\hfill \\ & {\displaystyle \underset{\mathsf{\Phi }}{\int }\mathsf{\delta }{\theta }_{\mathrm{d}}{K}_{\mathrm{R}}{\theta }_{\mathrm{d}}\rho (x)\mathrm{d}\varphi }-{\displaystyle \sum \mathsf{\delta }{\mathit{W}}^{\mathrm{T}}\mathit{Q}}-{\displaystyle \underset{\mathsf{\Phi }}{\int }\mathsf{\delta }{\mathit{W}}^{\mathrm{T}}\mathit{q}\rho (x)\mathrm{d}\varphi }=0\hfill \\ & \forall \mathsf{\delta }{\mathit{\epsilon }}_{\mathbf{s}},\mathsf{\delta }{\mathit{\epsilon }}_{\mathbf{c}},\mathsf{\delta }{\mathit{d}}_{\mathbf{slip}},\mathsf{\delta }{\theta }_{\mathrm{d}},\mathsf{\delta }\mathit{W}\hfill \end{array}$$

(9)

where ${\displaystyle \underset{\mathsf{\Phi }}{\int }{\displaystyle \underset{A_{\mathrm{s}}}{∯}\mathsf{\delta }{\mathit{\epsilon }}_{\mathbf{s}}^{\mathrm{T}}{\mathit{\sigma }}_{\mathbf{s}}\mathrm{d}a}\rho (x)\mathrm{d}\varphi }$ and ${\displaystyle \underset{\mathsf{\Phi }}{\int }{\displaystyle \underset{A_{\mathrm{c}}}{∯}\mathsf{\delta }{\mathit{\epsilon }}_{\mathbf{c}}^{\mathrm{T}}{\mathit{\sigma }}_{\mathbf{c}}\mathrm{d}a}\rho (x)\mathrm{d}\varphi }$ are the internal virtual work produced by the deformation of the steel beam and concrete slab, respectively; As and Ac denote the cross-sectional areas of the steel beam and concrete slab, respectively; and Φ is the central angle of the curved composite box beam. In terms of Equations (6) and (7), the strain variables of the steel beam and concrete slab can be written as follows:

$${\mathit{\epsilon }}_{\mathit{k}}=\mathit{S}{\mathit{B}}_{\mathit{k}}\mathit{d}$$

(10)

where the dimensions of the matrix B_k are 11 × 23, and its exact value is given in Appendix B. The matrix S and displacement vector d are expressed as follows:

$$\mathit{S}=\left[\begin{array}{ccccccccccc}\frac{R_0}{\rho (x)}& \frac{R_0}{\rho (x)}x& \frac{R_0}{\rho (x)}y& \frac{R_0}{\rho (x)}\omega & \frac{R_0}{\rho (x)}{\omega }_{\mathrm{d}}& \frac{R_0}{\rho (x)}{\psi }_{\mathrm{c}}& \frac{R_0}{\rho (x)}{\psi }_{\mathrm{s}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& \frac{R_0}{\rho (x)}{r}^{*}& \frac{\partial \omega }{\partial s}& {\psi }_{\mathrm{c},x}& {\psi }_{\mathrm{s},x}\end{array}\right]$$

(11)

$$\mathit{d}={\left\{\begin{array}{cccc}[\mathit{u}]& [\mathit{v}]& [\mathit{w}]& \begin{array}{cccc}[\mathit{\theta }]& [\mathit{\beta }]& [{\mathit{\theta }}_d]& \begin{array}{cc}[\mathsf{\Omega }]& [\mathit{f}]\end{array}\end{array}\end{array}\right\}}^{\mathrm{T}}$$

(12)

where

$$[\mathit{u}]=\left(\begin{array}{ccc}u& u^{\prime }& u^{″}\end{array}\right)$$

(13)

$$[\mathit{v}]=\left(\begin{array}{ccc}v& v^{\prime }& v^{″}\end{array}\right)$$

(14)

$$[\mathit{w}]=\left(\begin{array}{cc}w& w^{\prime }\end{array}\right)$$

(15)

$$[\mathit{\theta }]=\left(\begin{array}{cc}\theta & {\theta }^{\prime }\end{array}\right)$$

(16)

$$[\mathit{\beta }]=\left(\begin{array}{cc}\beta & {\beta }^{\prime }\end{array}\right)$$

(17)

$$[{\mathit{\theta }}_{\mathbf{d}}]=\left(\begin{array}{ccc}{\theta }_{\mathrm{d}}& {\theta }_{\mathrm{d}}{}^{\prime }& {\theta }_{\mathrm{d}}{}^{″}\end{array}\right)$$

(18)

$$[\mathsf{\Omega }]=\left(\begin{array}{cccc}2{\mathsf{\Omega }}_x& 2{\mathsf{\Omega }}_x{}^{\prime }& 2{\mathsf{\Omega }}_s& 2{\mathsf{\Omega }}_s{}^{\prime }\end{array}\right)$$

(19)

$$[\mathit{f}]=\left(\begin{array}{cccc}{f}_{\mathrm{c}}& f_{\mathrm{c}}{}^{\prime }& f_{\mathrm{s}}& {f^{\prime }}_{\mathrm{s}}\end{array}\right)$$

(20)

When the bridge is in use, the steel beam is in the elastic phase, and the stress–strain relationship of the steel beam is expressed as follows:

$${\mathit{\sigma }}_s=\left\{\begin{array}{c}{\sigma }_s\\ {\tau }_s\end{array}\right\}=E_s\left[\begin{array}{cc}1& 0\\ 0& 1/2\left(1+{\mu }_{\mathrm{s}}\right)\end{array}\right]\left\{\begin{array}{c}{\epsilon }_s\\ {\gamma }_s\end{array}\right\}$$

(21)

where σ_s and τ_s are the normal stress and shear stress of the steel beam, respectively; E_s is the elastic modulus of steel; and μ_s is the Poisson’s ratio of steel.

To take the effects of concrete shrinkage and creep into consideration, the stress–strain relationship of the concrete slab can be expressed as follows:

$$\left\{\begin{array}{c}{\epsilon }_c\left(t\right)-{\epsilon }_{\mathrm{c}}^{\mathrm{sh}}\left(t\right)\\ {\gamma }_c\left(t\right)\end{array}\right\}=J\left(t,t_0\right)\left[\begin{array}{cc}1& 0\\ 0& 2\left(1+{\mu }_{\mathrm{c}}\right)\end{array}\right]\left\{\begin{array}{c}{\sigma }_c\left(t_0\right)\\ {\tau }_c\left(t_0\right)\end{array}\right\}+{\displaystyle {\int }_{t_0}^{t}J\left(t,\tau \right)\left[\begin{array}{cc}1& 0\\ 0& 2\left(1+{\mu }_{\mathrm{c}}\right)\end{array}\right]\left\{\begin{array}{c}\mathrm{d}{\sigma }_c\left(\tau \right)\\ \mathrm{d}{\tau }_c\left(\tau \right)\end{array}\right\}}$$

(22)

To solve the integral problem in Equation (22), Bazant [21] proposed the age-adjusted effective modulus method to calculate the stress–strain relationship of concrete, which is expressed as follows:

$${\mathit{\sigma }}_{\mathbf{c}}\left(t\right)=\left\{\begin{array}{c}{\sigma }_{\mathrm{c}}\left(t\right)\\ {\tau }_{\mathrm{c}}\left(t\right)\end{array}\right\}=E_{\mathrm{ce}}\left(t,t_0\right)\left[\begin{array}{cc}1& 0\\ 0& 1/\left(2+2{\mu }_{\mathrm{c}}\right)\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{\mathrm{c}}\left(t\right)-{\epsilon }_{\mathrm{c}}^{\mathrm{sh}}\left(t\right)\\ {\gamma }_{\mathrm{c}}\left(t\right)\end{array}\right\}+\xi \left(t,t_0\right)\left\{\begin{array}{c}{\sigma }_{\mathrm{c}}\left(t_0\right)\\ {\tau }_{\mathrm{c}}\left(t_0\right)\end{array}\right\}$$

(23)

where t₀ is the initial loading age; ${\epsilon }_{\mathrm{c}}^{\mathrm{sh}}$ is the concrete shrinkage strain; μ_c is the Poisson’s ratio of concrete; J(t,t₀) is the concrete creep function at time t; σ_c and τ_c are the normal stress and shear stress of the concrete slab, respectively; E_ce(t,t₀) is the concrete age-adjusted effective elastic modulus; and ξ(t,t₀) is the influence factor of the creep effects. E_ce(t,t₀) and ξ(t,t₀) can be calculated by Equations (24) and (25), respectively.

$$E_{\mathrm{ce}}\left(t,t_0\right)=\frac{E_{\mathrm{c}}\left(t_0\right)}{1+\chi \left(t,t_0\right)\phi \left(t,t_0\right)}$$

(24)

$$\xi \left(t,t_0\right)=\frac{\phi \left(t,t_0\right)\left(\chi \left(t,t_0\right)-1\right)}{1+\chi \left(t,t_0\right)\phi \left(t,t_0\right)}$$

(25)

where E_c(t₀) is the concrete elastic modulus at time t₀, φ(t,t₀) is the concrete creep coefficient, and χ(t,t₀) is the concrete aging coefficient, which can be derived from Equation (26).

$$\chi \left(t,t_0\right)=\frac{E_{\mathrm{c}}\left(t_0\right)}{E_{\mathrm{c}}\left(t_0\right)-R\left(t,t_0\right)}-\frac{1}{\phi \left(t,t_0\right)}$$

(26)

where R(t,t₀) is the concrete relaxation function, which can be calculated according to the literature [21].

Note that ${\displaystyle \underset{\mathsf{\Phi }}{\int }{\displaystyle \underset{2b_{\mathrm{ts}}}{\int }\mathsf{\delta }{\mathit{d}}_{\mathbf{slip}}^{\mathrm{T}}{\mathit{q}}_{\mathbf{sh}}\mathrm{d}x}\rho (x)\mathrm{d}\varphi }$ in Equation (9) is the internal virtual work produced by interfacial slip between the steel beam and concrete slab. According to Equations (4) and (5), the interfacial slip between the steel beam and concrete slab can be written as follows:

$${\mathit{d}}_{\mathbf{slip}}={\left\{\begin{array}{cc}{u}_{\mathrm{sp}}& w_{\mathrm{sp}}\end{array}\right\}}^{\mathrm{T}}={\mathit{S}}_{\mathsf{\Omega }}{\mathit{B}}_{\mathsf{\Omega }}\mathit{d}$$

(27)

where the dimensions of matrix B_Ω are 2 × 23, as shown in Appendix B, and S_Ω is an element matrix with dimensions of 2 × 2.

During use, the shear stress flow of a unit area at the interface between the steel beam and concrete slab can be expressed as follows:

$${\mathit{q}}_{\mathbf{sh}}={\left\{\begin{array}{cc}{q}_{\mathrm{us}}& q_{\mathrm{ws}}\end{array}\right\}}^{\mathrm{T}}=\left[\begin{array}{cc}{\rho }_{\mathrm{u}}& 0\\ 0& {\rho }_{\mathrm{w}}\end{array}\right]\left\{\begin{array}{c}{u}_{\mathrm{sp}}\\ w_{\mathrm{sp}}\end{array}\right\}={\mathit{\rho }}_{\mathbf{slip}}{\mathit{d}}_{\mathbf{slip}}$$

(28)

where q_us and q_ws are transverse and longitudinal shear stress flow of a unit area in the composite beam interface, respectively; ρ_slip is the shear connection stiffness matrix of a unit area in the composite beam interface; and ρ_u and ρ_w are the transverse and longitudinal shear connection stiffness of a unit area in the composite beam, respectively.

Note that ${\displaystyle \underset{\mathsf{\Phi }}{\int }\mathsf{\delta }{\theta }_{\mathrm{d}}{K}_{\mathrm{R}}{\theta }_{\mathrm{d}}\rho (x)\mathrm{d}\varphi }$ in Equation (9) is the internal virtual work produced when the frame resists distortion. K_R is the frame distortion-resisting stiffness, which can be calculated according to the literature [22].

Note that ${\displaystyle \sum \mathsf{\delta }{\mathit{W}}^{\mathrm{T}}\mathit{Q}}$ and ${\displaystyle \underset{\mathsf{\Phi }}{\int }\mathsf{\delta }{\mathit{W}}^{\mathrm{T}}\mathit{q}\rho (x)\mathrm{d}\varphi }$ in Equation (9) are the external virtual work produced by external loads, wherein $\mathit{Q}={\left(\begin{array}{ccc}{Q}_n& Q_s& Q_z\end{array}\right)}^{\mathrm{T}}$ and $\mathit{q}={\left(\begin{array}{ccc}{q}_n& q_s& q_z\end{array}\right)}^{\mathrm{T}}$ are the concentrated load vector and the uniformly distributed load vector, respectively; Q_n, Q_s, and Q_z are concentrated forces in the directions of the n-axis, s-axis, and z-axis, respectively; q_n, q_s, and q_z are uniformly distributed forces in the directions of the n-axis, s-axis, and z-axis, respectively; and W is the displacement vector under external loads. The equation for this displacement vector is given in Equation (29), which is based on Equations (1)–(3).

$$\mathit{W}={\left\{\begin{array}{ccc}{\mathit{W}}_{nk}& {\mathit{W}}_{sk}& {\mathit{W}}_{zk}\end{array}\right\}}^{\mathrm{T}}={\mathit{H}}_1\mathit{A}+{\mathit{H}}_2{\mathit{A}}^{\prime }$$

(29)

where

$${\left[{\mathit{H}}_1\right]}_{3\times 10}={\left\{\begin{array}{ccc}{[{\mathit{B}}_1]}^{\mathrm{T}}& {[{\mathit{B}}_2]}^{\mathrm{T}}& {[{\mathit{B}}_3]}^{\mathrm{T}}\end{array}\right\}}^{\mathrm{T}}$$

(30)

$${\left[{\mathit{H}}_2\right]}_{3\times 10}={\left\{\begin{array}{ccc}{[\mathbf{0}]}_{1\times 10}{}^{\mathrm{T}}& {[\mathbf{0}]}_{1\times 10}{}^{\mathrm{T}}& {[{\mathit{B}}_4]}^{\mathrm{T}}\end{array}\right\}}^{\mathrm{T}}$$

(31)

$${\mathit{B}}_1=\left(\begin{array}{cccccccccc}\cos \alpha & -\sin \alpha & 0& -{\rho }_n& 0& -D_n& \frac{a_k\cos \alpha }{2}& 0& 0& 0\end{array}\right)$$

(32)

$${\mathit{B}}_2=\left(\begin{array}{cccccccccc}-\sin \alpha & -\cos \alpha & 0& {\rho }_s& 0& D_s& -\frac{a_k\sin \alpha }{2}& 0& 0& 0\end{array}\right)$$

(33)

$${\mathit{B}}_3=\left(\begin{array}{cccccccccc}0& 0& 1-x{k}_0& 0& -\omega & 0& 0& \frac{a_k\left(1-x{k}_0\right)}{2}& {\psi }_{\mathrm{c}}& {\psi }_{\mathrm{s}}\end{array}\right)$$

(34)

$${\mathit{B}}_4=\left(\begin{array}{cccccccccc}-x& -(y+\omega k_0)& 0& 0& 0& -{\omega }_{\mathrm{d}}& 0& 0& 0& 0\end{array}\right)$$

(35)

$$\mathit{A}={\left(\begin{array}{cccc}\begin{array}{cccc}u& v& w& \theta \end{array}& \beta & {\theta }_{\mathrm{d}}& \begin{array}{cc}2{\mathsf{\Omega }}_x& 2{\mathsf{\Omega }}_z{f}_{\mathrm{c}}{f}_{\mathrm{s}}\end{array}\end{array}\right)}^{\mathrm{T}}$$

(36)

2.5. Finite Beam Element for Curved Composite Box Beams

According to the methodology of finite element discretization, Equation (9) can be solved by adopting a two-node finite beam element, wherein thirteen DOFs exist at each of the nodes. The nodal displacement vector of the two-node finite beam element with 26 DOFs is expressed as follows:

$${\mathit{d}}^{\mathrm{e}}={\left\{\begin{array}{cc}{\mathit{d}}_1^{\mathbf{e}}& {\mathit{d}}_2^{\mathbf{e}}\end{array}\right\}}^T$$

(37)

$${\mathit{d}}_i^{\mathrm{e}}=\left(\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cc}{u}_i& u_i{}^{\prime }\end{array}& v_i& v_i{}^{\prime }& w_i\end{array}& {\theta }_i& {\beta }_i& {\theta }_{\mathrm{d}i}\end{array}& {\theta }_{\mathrm{d}i}{}^{\prime }& 2{\mathsf{\Omega }}_{xi}& \begin{array}{ccc}2{\mathsf{\Omega }}_{zi}& f_{\mathrm{c}i}& f_{\mathrm{s}i}\end{array}\end{array}\right)\mathrm{for}i=1,2$$

(38)

The shape function matrix [N]_{23 × 26} is introduced into Equation (12), and the shape function matrix [N_F]_{10 × 26} is introduced into Equation (36); both of the shape function matrices are shown in Appendix C. Substituting Equations (10)–(36) and the shape function matrices N and N_F into Equation (9), the finite beam element equation for the curved composite box beam can be derived as follows:

$$\mathit{K}{\mathit{d}}^{\mathrm{e}}=\mathit{F}$$

(39)

[K]_{26 × 26} in Equation (39) is the element stiffness matrix, which can be calculated as follows:

$$\mathit{K}={\displaystyle \underset{l_{\mathrm{e}}}{\int }{\mathit{N}}^{\mathrm{T}}\left[{\displaystyle \underset{A_s}{\iint }{\mathit{B}}_{\mathbf{s}}^{\mathrm{T}}{\mathit{S}}^{\mathrm{T}}{\mathit{E}}_{\mathbf{s}}\mathit{S}{\mathit{B}}_{\mathbf{s}}\mathrm{d}a}+{\displaystyle \underset{A_{\mathrm{c}}}{\iint }{\mathit{B}}_{\mathbf{c}}^{\mathrm{T}}{\mathit{S}}^{\mathrm{T}}{\mathit{E}}_{\mathbf{c}}\mathit{S}{\mathit{B}}_{\mathbf{c}}\mathrm{d}a}+{\displaystyle \underset{2b_{\mathrm{ts}}}{\int }{\mathit{B}}_{\mathsf{\Omega }}^{\mathrm{T}}{\mathit{S}}_{\mathsf{\Omega }}^{\mathrm{T}}{\mathit{\rho }}_{\mathbf{slip}}{\mathit{S}}_{\mathsf{\Omega }}{\mathit{B}}_{\mathsf{\Omega }}\mathrm{d}x+{\mathit{T}}_{\mathbf{R}}}\right]\mathit{N}\mathrm{d}z}$$

(40)

where [T_R]_{23 × 23} is the distortion-resisting stiffness matrix of the frame, in which T_R(13,13) = K_R and the other elements are 0, and l_e is the length of the finite beam elements.

[F]_{26 × 1} in Equation (39) is the equivalent load vector of the node, which can be calculated as follows:

$$\mathit{F}={\mathit{F}}_{\mathbf{q}}+{\mathit{F}}_{\mathsf{\xi }}+{\mathit{F}}_{\mathbf{sh}}$$

(41)

where F_q is the equivalent load vector of the node caused by external loads, which can be calculated by Equation (42); F_ξ is the equivalent load vector of the node caused by the concrete creep effect, which can be calculated by Equation (43); and F_sh is the equivalent load vector of the node caused by the concrete shrinkage effect, which can be calculated by Equation (44).

$${\mathit{F}}_{\mathbf{q}}={\displaystyle \underset{l_{\mathrm{e}}}{\int }\left({\mathit{N}}_{\mathbf{F}}^{\mathrm{T}}{\mathit{H}}_{\mathbf{1}}^{\mathrm{T}}+{{\mathit{N}}^{\prime }}_{\mathbf{F}}^{\mathrm{T}}{\mathit{H}}_{\mathbf{2}}^{\mathrm{T}}\right)\mathit{q}}\mathrm{d}z+\left({\mathit{N}}_{\mathbf{F}}^{\mathrm{T}}{\mathit{H}}_1^{\mathrm{T}}+{{\mathit{N}}^{\prime }}_{\mathbf{F}}^{\mathrm{T}}{\mathit{H}}_{\mathbf{2}}^{\mathrm{T}}\right)\mathit{Q}$$

(42)

$${\mathit{F}}_{\mathsf{\xi }}=-{\displaystyle \underset{l_{\mathrm{e}}}{\int }{\mathit{N}}^T{\mathit{B}}_{\mathbf{c}}{}^{\mathrm{T}}\left({\displaystyle \underset{A_{\mathrm{c}}}{∯}{\mathit{S}}^{\mathrm{T}}\xi \left(t,t_0\right){\mathit{\sigma }}_{\mathbf{c}}\left(t_0\right)\mathrm{d}a}\right)\mathrm{d}z}$$

(43)

$${\mathit{F}}_{\mathbf{sh}}={\displaystyle \underset{l_{\mathrm{e}}}{\int }{\mathit{N}}^T{\mathit{B}}_{\mathbf{c}}{}^{\mathrm{T}}\left({\displaystyle \underset{A_{\mathrm{c}}}{∯}{\mathit{S}}^{\mathrm{T}}{\mathit{E}}_{\mathbf{ce}}\left(t,t_0\right){\mathit{\epsilon }}_{\mathbf{c}}^{\mathbf{sh}}\left(t\right)\mathrm{d}a}\right)\mathrm{d}z}$$

(44)

where

$${\mathit{\sigma }}_{\mathbf{c}}\left(t_0\right)={\left\{\begin{array}{cc}{\sigma }_{\mathrm{c}}(t_0)& {\tau }_{\mathrm{c}}(t_0)\end{array}\right\}}^{\mathrm{T}}$$

(45)

$${\mathit{\epsilon }}_{\mathbf{c}}^{\mathbf{sh}}\left(t\right)={\left\{\begin{array}{cc}{\epsilon }_{\mathrm{c}}^{\mathrm{sh}}\left(t\right)& 0\end{array}\right\}}^{\mathrm{T}}$$

(46)

Notably, the calculation methods for the characteristic geometrical parameters of the section under torsion and distortion in the finite beam element model can be found in the literature [6,23].

3. Numerical Validation of the Beam Element Model

The proposed finite beam element model is applied to calculate the long-term behavior of a typical curved composite box beam. The cross-section and load layout of the curved composite box beam are shown in Figure 2: The arc length of the longitudinal section centroid line is 6200 mm, the central angle is 45°, and seven transverse diaphragms are uniformly distributed along the longitudinal direction; the slab is composed of C30 concrete and is 750 mm wide and 50 mm thick; the steel beam has a total height of 300 mm; the web has a thickness of 12 mm; the bottom flange is 350 mm wide and 12 mm thick; the interfacial shear connection stiffness along the transverse and longitudinal directions, ρ_u and ρ_w, are both equal to 100 N/mm³; and a vertical load of 10 kN/mm acts at the section centroid line.

The boundary conditions of the beam element model are as follows: At the mid-span, the longitudinal displacement and the longitudinal slip between the concrete slab and steel beam are all restrained to ensure symmetrical structural boundary constraints; at both ends of the beam, the transverse displacement, transverse slip, deflection, and torsion angle are restrained; and at the position of every diaphragm, the distortion angle is restrained.

The elaborate finite element model of a curved steel–concrete composite box beam is established in ANSYS. The concrete slabs, steel beams, and transverse diaphragms are modeled with Shell181 elements, whereas the shear connectors are modeled with Combin14 elements. The finite element model is shown in Figure 3.

The parameters of the material and environment of the proposed beam element model and the elaborate finite element model established in ANSYS are as follows: The elastic modulus of steel E_s and Poisson’s ratio of steel μ_s are 2.06 × 10⁵ MPa and 0.3, respectively; the calculation of the creep behavior of concrete in the two models use the age-adjusted effective modulus method; and the Poisson’s ratio of concrete μ_c is 0.2. The elastic modulus formula and the shrinkage creep model are obtained within CEB-FIP90 [24]; the drying age of concrete is three days, the relative humidity in the environment is 75%, and the initial loading age t₀ is 28 days. Using the finite beam element model and the elaborate finite element model, time-dependent changes in the deflection, transverse displacement, and stress of the curved composite box beam were calculated from the loading age to 3000 days, as shown in Figure 4, Figure 5 and Figure 6. From these three figures, the accuracy and applicability of the proposed two-node finite element with 26 DOFs are verified due to the good correlation between the results of the elaborate finite element model and those of the finite beam element model. Notably, the proposed beam element model and the elaborate FE model took 7 and 53 s to obtain the long-term behavior of the curved composite beam, respectively. It can be known that the calculation efficiency of the proposed beam element model is much higher than that of the elaborate FE model.

According to the three figures, the deflection and transverse displacement of the curved composite box beam increase over time due to concrete shrinkage and creep, and the compressive stress in the concrete decreases over time, whereas the tensile stress in the steel beam increases over time. Within 250 days after loading, the displacement and the stress of the curved composite box beam change significantly with respect to time, and then the growth rate slows and becomes gradual after the first year. The results show that the influence of concrete shrinkage and creep is significant in the early loading stage, whereas this influence is very small after the first year.

4. Parametric Analysis

4.1. Influence of the Time-Dependent Effects on Shear Lag

The shear-lag coefficient, λ, is introduced to reflect the influence of the shear-lag effect of curved composite beams. The shear-lag coefficient, λ, is defined as ratio of the stress at the intersection of the web and flange to that at the same location when the shear-lag effect is ignored. Based on the curved composite box girder in Section 3, which is used to verify the accuracy of the proposed model, the changes in the shear-lag coefficient of the concrete slab and steel bottom flange with respect to time are calculated by using the proposed finite beam element model, as shown in Figure 7 and Figure 8, respectively. The results show that the shear-lag coefficient at the inside and the outside of the concrete slab gradually diminishes over time. A small change in the shear-lag coefficient at the inside and outside of the steel bottom flange occurs in the curved composite beams over time, and it tends to stabilize in the late stage of development. The effects of concrete shrinkage and creep can weaken the shear-lag effects in curved composite box girders.

4.2. Influence of Time-Dependent Effects on Mechanical Behavior for Different Initial Curvatures

The initial central angles are set as 0°, 30°, and 60°, based on the curved composite box girder in Section 3, which is applied to verify the accuracy of the model. According to the proposed finite beam element model, the deflection, transverse displacement, and longitudinal displacement distributions along the span are calculated for a curved composite box girder at the initial and ultimate loading ages, as shown in Figure 9. The normal stress distributions along the transverse direction are calculated for the concrete slab and steel bottom flange at the initial and ultimate loading ages, as shown in Figure 10. Figure 9 shows that the maximum growth rates of deflection for central angles of 0°, 30°, and 60° are 30.5%, 29.2%, and 25.8%, respectively. The maximum growth rates of longitudinal displacement for central angles of 0°, 30°, and 60° are 20.3%, 20.6%, and 21.6%, respectively. The maximum growth rates of transverse displacement of the curved composite box girder with central angles of 30° and 60° are 47.2% and 45.1%, respectively. Figure 10 illustrates the stress distributions of the curved composite box girder with central angles of 0°, 30°, and 60° due to the influence of concrete shrinkage and creep. The tensile stress growth rates at the inside of the steel bottom flange for central angles of 0°, 30°, and 60° are 50.6%, 39.0%, and 18.6%, respectively. The compressive stress reduction rates at the inside of the concrete slab for central angles of 0°, 30°, and 60° are 58.8%, 57.8%, and 53.0%, respectively.

It can be concluded that the growth rate of displacement, caused by concrete shrinkage and creep effect, becomes small with the increase in the initial curvature of the curved box girder. The tensile stress growth rate of the steel girder decreases gradually with increasing initial curvature. The compressive stress reduction rate of the concrete slab decreases gradually with increasing initial curvature.

4.3. Influence of the Time-Dependent Effects on Mechanical Behavior for Different Interfacial Shear Connection Stiffness Values

The curved composite box girder in Section 3 was selected as the research object. The interfacial shear connection stiffness values are set as 0.1, 1, and 10 N/mm³. According to the proposed finite beam element model, the deflection and transverse and longitudinal displacement distributions along the span of the curved composite box girder at the initial and ultimate loading ages are calculated, as shown in Figure 11. The normal stress distributions along the transverse direction are calculated for the concrete slab and steel bottom flange at the initial loading age and ultimate loading age, as shown in Figure 12. Figure 11 shows that the maximum deflection growth rates are 54.4%, 32.3%, and 28.9% for curved composite box girders with shear connection stiffness values of 0.1, 1, and 10 N/mm³, respectively, due to the influence of concrete shrinkage and creep. Moreover, for shear connection stiffness values of 0.1, 1, and 10 N/mm³, the maximum growth rates of transverse displacement are 116.1%, 67.2%, and 49.2%, respectively, whereas the maximum growth rates of longitudinal displacement are 86.9%, 75.8%, and 71.3%, respectively. Figure 12 shows that the tensile stress growth rates at the inside of the steel bottom flange are 47.1%, 35.2%, and 37.8% for curved composite box girders with shear connection stiffness values of 0.1, 1, and 10 N/mm³, respectively, due to the influence of concrete shrinkage and creep. Under these same stiffness values, the compressive stress reduction rates at the inside of the concrete slab are 44.6%, 59.4%, and 57.1%, respectively.

It can be concluded that the growth rate of displacement, which is caused by concrete shrinkage and creep, becomes small as the interfacial shear connection stiffness of the curved composite box girder increases. The tensile stress growth rate of the steel girder decreases with increasing interfacial shear connection stiffness. The compressive stress reduction rate of the concrete slab becomes large with increasing interfacial shear connection stiffness.

5. Conclusions

(1) Based on the Vlasov’s one-dimensional theoretical model of curved beams, a torsional warping function, a distortional angle, transverse and longitudinal slip functions, a shear-lag intensity function for the concrete slab, a shear-lag intensity function for the steel beam, and constitutive relationships for concrete shrinkage and creep were introduced. The variation of the curvature in the radial direction was also taken into account. The equilibrium equation for the curved composite box beam was proposed according to the virtual work principle. Then, the stiffness matrix and equivalent nodal load vector of the two-node finite beam element with 26 DOFs, which analyzes the long-term behavior of curved composite box beams, were obtained by adopting the finite element method and the age-adjusted effective modulus method.

(2) An elaborate finite element model of a typical curved composite box beam was established. The correctness and applicability of the two-node finite beam element with 26 DOFs were verified through a comparison of the results from the finite beam element model and those from the elaborate finite element model.

(3) The proposed finite beam element model was adopted to analyze the influence of the concrete time-dependent effects on the long-term behavior of curved composite box beams. The analysis shows that significant changes of the displacement, stress, and shear-lag coefficient occur in the curved composite beams within the first year of loading, due to the time-dependent effects of concrete, after which the variation tendency becomes gradual.

(4) The proposed finite beam element model was adopted to analyze the influence of the concrete shrinkage and creep effects on the long-term behavior of curved composite box beams under different central angles and shear connection stiffness. The results show that increases of the central angle and the shear connection stiffness reduce the change rates of both displacement and stress, with respect to time.