A Beam Finite Element Model Considering the Slip, Shear Lag, and Time-Dependent Effects of Steel–Concrete Composite Box Beams

: A beam ﬁnite element model considering the slip, shear lag, and time-dependent effects of steel–concrete composite box beams have been proposed in this study. The element is employed to a one-dimensional analytical method that is solved involving an expression of spatial and time-dependent variables. A step-by-step method that does not involve storing the stress and strain histories, which is more accurate than the single-step algebraic method, is employed to solve the time variables. A recursive method was elaborated to determine spatial and time-dependent variables through the above method. The validity of the proposed method in instantaneous analysis is attested by the numerical data of the elaborate ﬁnite element model established in a commercial software, ANSYS, and that of the time-dependent analysis is veriﬁed by the existing test results on the long-term performance of composite beams. The proposed beam ﬁnite element model is applied to predict the time-related behavior of simply supported composite beams after validation. The results show that concrete shrinkage and creep signiﬁcantly inﬂuence the structural responses of the composite box beams. From the initial load on the 28th day to that in the 3rd year, the vertical deﬂection at the cross-section of the mid-span increased by 47.01%. The interface slip at the end increased by − 10.99%. The warping intensity function of the concrete slab and the steel beam at the end caused by shear lag increased by 111.64% and 7.01%, respectively. The maximum compressive stress on the concrete slab and the maximum tensile stress at the steel bottom ﬂange increased by − 6.75% and 4.56%, respectively.


Introduction
Steel-concrete composite beams, by far, have been employed into buildings and bridges world-wide. Design analysis and computational modeling of composite beams have long been active study areas [1][2][3][4][5]. Their concrete slab and steel beam are subjected to pressure and tension, respectively, under vertical loading, fully utilizing the good compressive strength of concrete and the excellent tensile strength of steel. The shear connectors serve as a kind of connection between the steel beam and the concrete slab to achieve the overall working performance. However, deformation of the shear connectors causes interface slip between the concrete slab and the steel beam, which reduces the bending stiffness of the composite beams. Actually, the interface slip concerns are characteristic not only of classical composite beams, but also of non-standard composite beams such as hybrid trussed beams with inclined shear studs [6,7]. Additionally, the shear lag effect on the slab causes varying stress distribution, remarkably for composite beams with an extremely wide concrete slab where the shear lag is positive or negative [8][9][10][11]. In the positive case, the stress at the slab-beam intersection is more than that for the remainder of the slab, while in the negative case, the stress is lower than that for the remainder of the slab. Due to the non-uniform distribution of stresses, disregarding the shear lag may bring out underestimation of the actual stresses in parts of the concrete slab. Subsequently, to simulate the mechanical behavior of composite beams accurately, the shear lag and interface slip must be taken into consideration. Figure 1 shows the two mechanical behaviors of the composite box beam above. On top of the mentioned spatial kinematic effects, the time-dependent effect of the structure because of the shrinkage and creep of the concrete slab have been a major focus in the study of composite beams, where creep and shrinkage of concrete lead the extra internal force and deformation of the composite beam to continuously change as time goes by [12,13]. Some studies have focused on the partial time-dependent effects of the composite beams. For example, D Huang et al. [14] studied the time-dependent shear lag of assembled composite beams with the concrete age difference across the section. Souici et al. [15] proposed a method to study the effects of creep in composite beams, which assumes a perfect connection between the concrete slab and the steel beam. Based on generalized beam theory, Henriques et al. [16] proposed a finite element to analyze the time-dependent effects of steel-concrete composite beams without considering the slip. face slip must be taken into consideration. Figure 1 shows of the composite box beam above. On top of the mentione time-dependent effect of the structure because of the shrin slab have been a major focus in the study of composite bea of concrete lead the extra internal force and deformation o uously change as time goes by [12,13]. Some studies hav dependent effects of the composite beams. For example, D time-dependent shear lag of assembled composite beams w across the section. Souici et al. [15] proposed a method t composite beams, which assumes a perfect connection bet steel beam. Based on generalized beam theory, Henriques ment to analyze the time-dependent effects of steel-conc considering the slip.  One-dimensional models that consider interface slip, behaviors [17][18][19][20][21] have been widely utilized in the behavio The one-dimensional beam model retains the longitudinal the cross-sectional characteristics are simplified in the tran which is obviously more efficient than the three-dimensi rectness and applicability of the one-dimensional model plenty of engineering practices.
The spatial and time-dependent variables are predict models. To solve the spatial variables, displacement functi of the entire beam are substituted into equations about v method (FEM) [20] is often utilized to solve the equations. studies on finite element models of GFRP composite beam variables is the main focus of the model. Many studies hav ent behavior of composite beams due to the creep and shr the shrinkage, since the shrinkage has little to do with the applied to the structure during the calculation, which is sim this is more complicated since the concrete creep is inext One-dimensional models that consider interface slip, shear lag, and time-dependent behaviors [17][18][19][20][21] have been widely utilized in the behavior of composite beam's analyses. The one-dimensional beam model retains the longitudinal direction (along the span) and the cross-sectional characteristics are simplified in the transverse and vertical directions, which is obviously more efficient than the three-dimensional elaborate model. The correctness and applicability of the one-dimensional model have been widely verified by plenty of engineering practices.
The spatial and time-dependent variables are predicted to analyze one-dimensional models. To solve the spatial variables, displacement functions and the intensity functions of the entire beam are substituted into equations about virtual work. The finite element method (FEM) [20] is often utilized to solve the equations. Madenci et al. conducted some studies on finite element models of GFRP composite beams [22][23][24][25]. How to consider time variables is the main focus of the model. Many studies have considered the time-dependent behavior of composite beams due to the creep and shrinkage of concrete [26][27][28][29][30][31][32][33][34][35][36]. For the shrinkage, since the shrinkage has little to do with the stress state, the initial strain is applied to the structure during the calculation, which is simpler to address. For the creep, this is more complicated since the concrete creep is inextricably relevant to the current stress state, and the stress state in the whole stress history, which explains why many studies focus on it.
Common solutions for time variables in the literature are general step-by-step methods and methods using single-step algebraic equations. The time-integration procedures are an essential part in numerical solution procedures that translate genetic integration relations into constitutive relations for the time decomposition that can be easily addressed in the solution algorithm. Methods using one-step algebraic equations, such as the effective modulus (EM) [26], mean stress (MS) [27], and age-adjusted effective modulus (AAEM) [28,29] methods, can be used for different varieties of orthogonal formulas used for numerical integration, disregarding stress histories, but they lose the accuracy. Instead, Bazant [30] proposed a general step-by-step method aimed at this numerical integration based on the trapezoidal or midpoint rule. The drawback, however, is that the stress history as a whole needs to be restored during the solving process, which leads to a long solving time and high computational cost. However, the problem has been subsequently solved by Bazant and Wu [31]. They proposed a new strategy for the step-by-step integration methods by introducing the Dirichlet series to fit the creep function, which only requires storing the stress and strain of the previous steps in time histories. This work paves a good numerical way for the accurate analysis of the structural time-dependent behavior while saving memory and improving computational efficiency.
The above findings show that without storing data on the stress history, the stepby-step method in the time history and the finite beam element method with additional degrees of freedom (DOFs) function based on the Euler Bernoulli beam in the space provide the most accurate prediction of the complex spatial time-dependent behavior of composite beams. The novelty of the study is to propose a one-dimension theoretical model of a composite box beam considering slip, shear lag, and time-dependence based on the virtual work principle. For the solution to this model, a finite element discretization method is employed for the spatial variables, and an incremental step-by-step method without storing stress and strain histories is utilized for the time variables. A beam element with 18 DOFs for the steel-concrete composite box beam is proposed. The accuracy of the proposed beam model is verified by numerical simulations and some classical tests. Among them, the instantaneous behavior of the model is verified by numerical calculations using the finite element model, and the time-dependent behavior of the model is verified by long-term performance test results.

Principal Knowledge
The beam element model was built on the basis of the classical Euler Bernoulli beam theory and the virtual work principle. Considering the behaviors of the steel-concrete composite box beam, the following basic assumptions are made for the theoretical derivation and modeling: (i) The vertical bending curvatures and transverse bending curvatures of the concrete slab and steel beam are identical; (ii) The deflections of the concrete slab and steel beam in both vertical and transverse deflections are identical; (iii) Shear deformation of the beam caused by bending is disregarded; (iv) The slip between the steel beam and the concrete slab is considered only in a longitudinal direction. The shear connections are arranged uniformly along the span so that the shear connection stiffness of the interface remains identical along the span; (v) The shear lag on vertical deflection is considered; (vi) The study only focused on the structure in the normal service stage, so the concrete slab is always under elastic stage. The concrete creep is simulated by use of a linear creep model according to Bazant's study. (vii) The study only focused on the structure in the normal service stage, so the steel beam is always under the elastic stage; (viii) The study only focused on the mechanical properties in the normal service stage, so the shear connections are always under elastic stage. Figure 2 shows the definition of the geometrical parameters of the I-beam or box combination beam, where O c is the centroid of the concrete slab, O s is the centroid of the steel beam, and C s is the torsional centroid of the transformed section of the composite beam.

Kinematics of the Composite Box Beam
From the above assumptions and the coordinate system in Figure 2, the transverse displacement u(x, y, z) and vertical displacement v(x, y, z) at any point of the beam, displacement of the concrete slab wc(x, y, z), and displacement of the steel beam ws(x, y, z) in longitudinal direction are expressed as follows: where yh is the y-directional coordinate of Cs; xh is the x-directional coordinate of Cs; wc0 and ws0 are the longitudinal displacements of Oc and Os, respectively; u0 is the transverse displacement of Oc or Os; v0 is the vertical displacement of Oc or Os; ϕis the torsion angle of the overall structure; fc and fs are the intensity functions of warping due to shear lag on the concrete slab and steel beam, respectively; and ψc(x) and ψs(x) are the functions of warping due to shear lag on the concrete slab and steel beam, respectively. The functions of warping on the concrete slab are calculated based on Equation (2). For the steel U-beam, there is no warpage due to shear lag in the upper flange and web; thus, ψs(x) = 0. For the bottom flange of the steel beam, the shear lag needs to be considered and is calculated by The magnitude of interface slippage between the steel beam and the concrete slab dsl is:

Kinematics of the Composite Box Beam
From the above assumptions and the coordinate system in Figure 2, the transverse displacement u(x, y, z) and vertical displacement v(x, y, z) at any point of the beam, displacement of the concrete slab w c (x, y, z), and displacement of the steel beam w s (x, y, z) in longitudinal direction are expressed as follows: where y h is the y-directional coordinate of C s ; x h is the x-directional coordinate of C s ; w c0 and w s0 are the longitudinal displacements of O c and O s , respectively; u 0 is the transverse displacement of O c or O s ; v 0 is the vertical displacement of O c or O s ; φ is the torsion angle of the overall structure; f c and f s are the intensity functions of warping due to shear lag on the concrete slab and steel beam, respectively; and ψ c (x) and ψ s (x) are the functions of warping due to shear lag on the concrete slab and steel beam, respectively. The functions of warping on the concrete slab are calculated based on Equation (2). For the steel U-beam, there is no warpage due to shear lag in the upper flange and web; thus, ψ s (x) = 0. For the bottom flange of the steel beam, the shear lag needs to be considered and is calculated by Equation (3).
The magnitude of interface slippage between the steel beam and the concrete slab d sl is: where h 0 is the vertical distance between O c and O s . From the variables of the displacement of the concrete slab, the normal strain ε c and shear strain γ c of the concrete slab are obtained: From the displacement variables of the steel beam, the normal strain ε s and tangential strain γ s of the steel beam are obtained: where r * c is the vertical distance from C s to any point of the concrete slab along its thin wall and r * s is the vertical distance from C s to any point of the steel beam along its thin wall of. ψ c,x is the derivative of ψ c of x in first order, and ψ s,x is the first-order derivative of ψ s with respect to x.

Virtual Work of the Composite Box Beam
The one-dimensional analysis model of the composite box beam is on the basis of the virtual work principle where the expression is: where A s and A c are the areas of the steel beam and concrete slab cross-section, respectively, and L is the length of the composite beam span. In order to present the results concisely, the variables and equations are subsequently expressed as a matrix.
(1) Internal virtual work of the steel beam and the concrete slab δε T s σ s dadz in Equation (7) represent the internal virtual work caused by the deformation of the concrete slab and the steel beam, separately. ε c is the strain matrix of the concrete slab, as shown in Equation (8), which includes the normal strain and tangential strain, respectively.
The strain matrix ε c of a concrete slab is: where ε s is the strain matrix of the steel beam, as shown in Equation (12), including normal strain and tangential strain, respectively.
The strain matrix of steel beam ε s is: where The stress matrix of the concrete slab σ c (Equation (16)) including normal stress σ c and tangential stress τ c : Based on assumption (vi), considering the shrinkage and creep effect of concrete, introduce the creep function J(t,t 0 ) and shrinkage strain matrix ε c,sh (t). The relationship between the stress and the strain of concrete is expressed as follows: where t is the age of concrete, t 0 is the initial loading age of concrete, and ε c,sh is the matrix of instantaneous shrinkage strain of concrete with respect to t (Equation (18)). ε c,sh (t) in Equation (18) is the instantaneous shrinkage strain of concrete with respect to t. The relationship X c is shown in Equation (19), and µ c in Equation (19) is Poisson's ratio of concrete. ε c,sh = ε c,sh 0 T (18) The stress matrix of the steel beam σ s (Equation (20)) comprises the normal stress σ s and shear stress τ s .
According to assumption (vii), the relationship between σ s and ε s is: where the relationship X s is shown in Equation (22), E s is the elastic modulus of steel, and µ s is Poisson's ratio of steel.
(2) Internal virtual work by interface slip L δd sl q sl dz in Equation (7) is the internal virtual work due to the interface slip between the steel beam and the concrete slab. From the assumption (viii), the relationship between the interface shear force d sl and the interface slip displacement q sl is expressed as follows: where ρ is the per unit length along the longitudinal direction for the shear connection stiffness of the interface.

(3) External virtual work
∑ δW T Q in Equation (7) is the external virtual work due to the concentrated load at any position of the composite beam, and L δW T qdz is the external virtual work due to the distributed load at any position of the composite beam.
According to Equation (1), the displacement vector W at any position of the steel beam or the concrete slab is expressed as: Furthermore, Equation (24) is written as: where The concentrated load Q and distributed load q are expressed as follows: where Q x , Q y , and Q z are the concentrated loads in the x, y, and z direction, respectively, and q x , q y , and q z are the distributed loads in the x, y, and z direction, respectively. It should be noted that this analytical model considers the case that the load acts at any position.

Numerical Procedure of the Analysis Model
The one-dimensional theoretical model of the composite box beam includes two parts as follows. First, the time-dependent components are obtained using a more accurate step-by-step method. To save computational space (without storing the stress history), the Dirichlet series is introduced to the creep function. Second, the space-related components are obtained by the FEM selected to discretize the composite box beam into several beam elements with two nodes. Next, the components of element stiffness and equivalent load matrices at each node are calculated, which are applied to obtain the overall stiffness matrix and the overall equivalent load matrix.

Time Integration: Incremental Step-by-Step Method without Storing the Histories of Stress and Strain
When the effect of temperature on concrete is disregarded, the gross strain of concrete is the combination of strain, shrinkage strain, and creep strain in an instant [34], as shown in Equation (30). The shrinkage strain is not related to the forces of the structure, and the numerical model of the overall structure equals to the initial strain of the composite box beam when considering shrinkage, which is a relatively simple method to address. The instantaneous strain and creep strain are related to the forces of the structure. Moreover, the strain is related to the stress, both of which change over time. So, the numerical analysis of this process is far more cumbersome. According to reference [34], the creep function J(t, t 0 ) is employed to obtain the instantaneous strain and creep strain.
According to reference [34], the creep function is written in the form of Equation (31): where E c (t) is the elastic modulus of concrete in each time step t, and C(t,t 0 ) is the specific creep function of concrete in each time step t when the initial loading age is t 0 . By substituting Equation (31) into Equation (17) and separating the instantaneous strain ε c,e (t) and creep strain ε c,cr (t), the following expressions are obtained: To improve the accuracy as much as possible, the incremental step-by-step method is adopted, and the Dirichlet series is employed into the specific creep function C(t, t 0 ) by the step-by-step method (without storing the stress and strain histories, respectively). According to reference [34], the Dirichlet series expansion of C(t, t 0 ) in Equation (31) is shown in Equation (34), i.e., the Kabir formula: Bazant [34] suggested that for the concrete creep problem, using m = 4 is sufficient to guarantee the validity of the expansion of series, and the delay time τ j = 10 j−1 (j = 1, 2, . . . , m). The α j (t 0 ) (i = 1, 2, . . . , m) in Equation (34) is solved by the least-squares method provided that C(t, t 0 ) is known.
To apply the incremental step-by-step method, the time step from t 0 to t is discretized in the form of ∆t n = t n − t n−1 (n ≥ 1), and ∆ε n c , ∆ε n c,cr , ∆ε n c,sh and ∆σ n c of ∆t n are solved. Among these variables, ∆ε n c,sh is independent of stress and is obtained from the timedependent constitutive relation of shrinkage strain; ∆ε n c is obtained from the displacement matrix d c of concrete, while the calculation of ∆ε n c,cr and ∆σ n c is more complicated and needs to be calculated by the recursive method [34]: Using the nth time step ∆t n as an example, we have Equation (33): When n > 0, When Therefore, In Equation (38), Therefore, we obtain: It can be determined that β n j satisfies the recurrence relation, as shown in Equation (42). It is only necessary to store β n−1 j and stress increment ∆σ n−1 c at the previous time step ∆t n−1 to obtain β n j at time step ∆t n , which avoids the need to store stress and strain histories.
Therefore, it is determined that the strain increment ∆ε n c,cr of time step ∆t n is obtained from the stress increment ∆σ n−1 c of the previous time step ∆t n−1 and β n j of the time step ∆t n . According to Equation (32), the following expressions are obtained: where As shown in Equation (43), the ∆σ n c of time step ∆t n is obtained. Afterward, one can proceed to the next time step ∆t n+1 and continue to repeat the above calculation process.

Space Integration: The Finite Beam Element with 18 DOFs
An FEM with high accuracy and stability is applied to obtain the one-dimensional analysis model in space. Like Equation (7), the identical virtual work principle is followed in incremental time steps ∆t n : The composite box beam is discretized into two-node finite beam elements with 18 DOFs, each node containing 9 DOFs. The elements of the node displacement matrix of d e are:  (45), and the equilibrium equation of the composite box beam at time interval ∆t n is solved as: where K n is the incremental stiffness matrix for the incremental time interval ∆t n , ∆d n e is the incremental node displacement matrix of the element for the time incremental interval ∆t n , and ∆F n is the incremental equivalent load matrix of the element node for the time incremental interval ∆t n .
The form of the expression for the incremental stiffness matrix K n is: where l e is the length of the finite beam element. The elements of incremental equivalent node load matrix of ∆F n are calculated by the following expression: ∆F n consists of three parts, where ∆F n ext is the equivalent node load matrix externally produced in the incremental step, ∆F n cr is the equivalent node load matrix of produced by the creep strain in the incremental step, and ∆F n sh is the equivalent node load matrix caused by the shrinkage strain in the incremental step. The equations of the three load matrices above are: Combining the above calculations in the time and space domains, a solution procedure for a two-node finite beam element model with 18 DOFs of the composite box beam considering slip, shear lag, and time-dependence is proposed: First, based on Equation (34) and C(t n ,t i ) (i = 0, 1, . . . , n−1) for all of the time steps, α j (t i ) are obtained by least-square fitting, after which the recursive process begins: (i) Calculate β n j according to Equation (42) (48) and (9); (v) Calculate ∆σ n c according to Equation (43) for the next incremental time step; (vi) Return to (i) and perform the iterative computation for the next incremental time step. Figure 3 shows the synoptic scheme of inputted geometric and mechanical parameters for the solution of the proposed beam element model. Figure 4 shows the flow of the solution procedure for the proposed beam element.    In short, this study proposed a finite beam element with 18 DOFs for a steel-concrete composite box beam considering interface slip, shear lag, and time-dependence.

Verification of the Beam Finite Element Model
To verify the accuracy and applicability of the proposed beam finite element model, the instantaneous and long-term results of the finite beam element model were verified.    In short, this study proposed a finite beam element with 18 DOFs for a steel-concrete composite box beam considering interface slip, shear lag, and time-dependence. In short, this study proposed a finite beam element with 18 DOFs for a steel-concrete composite box beam considering interface slip, shear lag, and time-dependence.

Verification of the Beam Finite Element Model
To verify the accuracy and applicability of the proposed beam finite element model, the instantaneous and long-term results of the finite beam element model were verified.

Instantaneous Behavior
First, the instantaneous behavior of a simply supported composite beam applied to a project under the action of a uniform load is analyzed. The standard span length of the composite beam is 60 m, and the parameters of the cross-section are illustrated in Table 1.  The elaborate finite element (FE) model and the proposed finite beam element model of the simply supported composite beam were utilized to numerically analyze the structure. The elaborate FE model was established by using a commercial software, ANSYS version 18.0, as follows: both the concrete slab and steel beam were meshed by using shell elements (four-node general order Shell181 element in ANSYS). The studs at the interface between the steel plate and the concrete were simulated by using a linear spring element (two-node general order Combin14 element in ANSYS). The elastic stiffness of the shear studs was calculated as per CEB-FIP Model Code 1990 [37]. Figure 5 shows the elaborate FE model of the simply supported composite beam in ANSYS, which has 4860 nodes and 5612 elements. The mesh size of the elaborate element model is approximately 600 mm.  The beam element was established by using MATLAB 2019. A mesh sensitivity test was performed in which 20 beam elements and 21 nodes were used to mesh the specimens. The results showed that the sensitivity of the FE simulation results is within 1% of the observed values. The boundary conditions for the simply supported beam are listed in Table 2. In order to ensure simulation of the boundary condition of the simply supported beam, the longitudinal displacements of concrete slab and steel girder wc and ws and the intensity functions of warping due to shear lag on concrete slab and steel girder fc and fs are restricted in the mid-span. The transverse displacement u, vertical displacement v, and torsion angle ϕ are restricted at both ends of the beam.

Location of Section
Restraint DOFs Both ends of the beam u, v, ϕ The beam element was established by using MATLAB 2019. A mesh sensitivity test was performed in which 20 beam elements and 21 nodes were used to mesh the specimens. The results showed that the sensitivity of the FE simulation results is within 1% of the observed values. The boundary conditions for the simply supported beam are listed in Table 2. In order to ensure simulation of the boundary condition of the simply supported beam, the longitudinal displacements of concrete slab and steel girder w c and w s and the intensity functions of warping due to shear lag on concrete slab and steel girder f c and f s are restricted in the mid-span. The transverse displacement u, vertical displacement v, and torsion angle φ are restricted at both ends of the beam.

Location of Section Restraint DOFs
Both ends of the beam u, v, φ Mid-span w c , w s , f c , f s Figure 6 shows the calculation results of the elaborate FE model and proposed finite beam element model, including vertical deflections, interface slip, stresses on the concrete slab, and stresses at the bottom flange of the steel beam. The figure shows that the deviation between the two calculations is only within 5%.
the observed values. The boundary conditions for the simply supported beam are listed in Table 2. In order to ensure simulation of the boundary condition of the simply supported beam, the longitudinal displacements of concrete slab and steel girder wc and ws and the intensity functions of warping due to shear lag on concrete slab and steel girder fc and fs are restricted in the mid-span. The transverse displacement u, vertical displacement v, and torsion angle ϕ are restricted at both ends of the beam.

Location of Section
Restraint DOFs Both ends of the beam u, v, ϕ Mid-span wc, ws, fc, fs Figure 6 shows the calculation results of the elaborate FE model and proposed finite beam element model, including vertical deflections, interface slip, stresses on the concrete slab, and stresses at the bottom flange of the steel beam. The figure shows that the deviation between the two calculations is only within 5%. In terms of the computation cost, the elaborate FE model requires 15.32 s, while the proposed beam element requires a total of 1.56 s. Therefore, the proposed beam element has notably higher computational efficiency than the elaborate FE model.

Time-Dependent Behavior
Fan et al. [38] carried out a test of two composite beams under vertical uniform loads in the long term. Both specimens have identical geometry configuration: the span of the composite box beam is 4 m, the steel beam height is 180 mm, the thickness of the reinforced concrete slab is 60 mm, and its width is 600 mm. A row of studs was welded on the upper flange of the steel beam with an interval of 80 mm longitudinally. The concrete slabs of In terms of the computation cost, the elaborate FE model requires 15.32 s, while the proposed beam element requires a total of 1.56 s. Therefore, the proposed beam element has notably higher computational efficiency than the elaborate FE model.

Time-Dependent Behavior
Fan et al. [38] carried out a test of two composite beams under vertical uniform loads in the long term. Both specimens have identical geometry configuration: the span of the composite box beam is 4 m, the steel beam height is 180 mm, the thickness of the reinforced concrete slab is 60 mm, and its width is 600 mm. A row of studs was welded on the upper flange of the steel beam with an interval of 80 mm longitudinally. The concrete slabs of specimens LCB1 and LCB2 were constructed of C20 concrete and C30 concrete, respectively. The cubic compressive strengths on the 7th day and the 28th day of the concrete applied in LCB1 were tested to be 24.3 MPa and 32.3 MPa, respectively. The cube compressive strengths on the 7th day and 28th day of the concrete applied in LCB2 were 33.4 MPa and 44.7 MPa, respectively. The composite beams were loaded from the 7th day of age continuously and constantly and had been monitored for 3 years. The vertical uniform load (including the weight of the beam itself) was 6.23 kN/m. Bradford and Gilbert [39] conducted long-term performance tests on four simply supported steel-concrete composite beams with specimens named B1 to B4. Specimens B1 and B2 arrange two rows of studs with an interval 200 mm longitudinally. Specimens B3 and B4 arrange two rows of studs with an interval 600 mm longitudinally. The uniform vertical loads are sustained by the specimens B2 and B4 considering self-weight (1.92 kN/m) only, while those that are sustained by the specimens B1 and B3 consider selfweight (1.92 kN/m) and extra uniform loads (7.52 kN/m). The cylindrical compressive strength and elastic modulus of the concrete applied in the specimens is 31.1 MPa and 2.51 × 10 4 MPa, respectively. The elastic modulus of steel is 2.0 × 10 5 MPa, the final value of shrinkage of concrete on the 220th day was 410 × 10 −6 , and the coefficient of creep was 2.6.  Bradford and Gilbert [39] conducted long-term performance tests on four simply supported steel-concrete composite beams with specimens named B1 to B4. Specimens B1 and B2 arrange two rows of studs with an interval 200 mm longitudinally. Specimens B3 and B4 arrange two rows of studs with an interval 600 mm longitudinally. The uniform vertical loads are sustained by the specimens B2 and B4 considering self-weight (1.92 kN/m) only, while those that are sustained by the specimens B1 and B3 consider self-weight (1.92 kN/m) and extra uniform loads (7.52 kN/m). The cylindrical compressive strength and elastic modulus of the concrete applied in the specimens is 31.1 MPa and 2.51 × 10 4 MPa, respectively. The elastic modulus of steel is 2.0 × 10 5 MPa, the final value of shrinkage of concrete on the 220th day was 410 × 10 −6 , and the coefficient of creep was 2.6. Figure 8 displays the measured vertical displacements of the specimens and the predicted results of the proposed finite beam element model. As for the model, the coefficient curves of the shrinkage strain and creep coefficient of concrete are determined by the material tests. The prediction results of the proposed model show agreement with the test results. × 10 MPa, respectively. The elastic modulus of steel is 2.0 × 10 MPa, the final value of shrinkage of concrete on the 220th day was 410 × 10 −6 , and the coefficient of creep was 2.6. Figure 8 displays the measured vertical displacements of the specimens and the predicted results of the proposed finite beam element model. As for the model, the coefficient curves of the shrinkage strain and creep coefficient of concrete are determined by the material tests. The prediction results of the proposed model show agreement with the test results. A comparison between the experimental data given in [38,39] and the prediction model proposed by Zhu and Su [35] was carried out. The comparison results obtained by Zhu and Su are similar to those obtained by this study, as shown in Figures 7 and 8. A comparison between the experimental data given in [38,39] and the prediction model proposed by Zhu and Su [35] was carried out. The comparison results obtained by Zhu and Su are similar to those obtained by this study, as shown in Figures 7 and 8.

Application of the Beam Finite Element Model
The proposed method is applied to predict the time-dependent behaviors of simply supported composite beams. The span of the composite beam is 40 m, the parameters of the cross-section shown in Table 3. The material properties of the specimens are as follows: shear connection stiffness ρ = 1 kN/mm 2 , elastic modulus of steel E s = 2.06 × 10 5 MPa, Poisson's ratio of steel µ s = 0.3, cubic compressive strength of concrete f ck = 50 MPa, elastic modulus of concrete E c = 3.8629 × 10 4 MPa, and Poisson's ratio of concrete µ c = 0.2, the shrinkage age t sh = 7 d, loading age t 0 = 28 d, and relative humidity RH = 75%. A vertical uniform load of 200 kN/m was applied to the composite beam. The shrinkage strain and creep coefficient calculation methods specified in CEB-FIP 90 [37] were adopted in the model.

Influence of the Shear Connection Stiffness
To investigate the influence of shear connection stiffness, two other composite beam models with dissimilar shear connection stiffnesses were also analyzed for the cases ρ = 0.5 kN/mm 2 and ρ = 10 kN/mm 2 . Figures 9-12 show the vertical deflections, interface slip, the stresses of the concrete slab, and the bottom flange of the steel beam, where (a) and (b) represent the 28th day and 3rd year, respectively. According to the above figures, compared to the model with ρ = 1 kN/mm 2 , the mid-span vertical deflections of the model with ρ = 0.5 kN/mm 2 on the 28th day and in the 3rd year are 13.90% larger and 8.40% larger, respectively; for the interface slip at the beam end, the vertical deflections are 91.64% larger and 90.62% larger, respectively; for the maximum compressive stress on the concrete slab, the vertical deflections are 2.2% smaller and 2.13% smaller, respectively; for the maximum tensile stress in the steel beam bottom flange, the vertical deflections are 1.06% larger and 0.96% larger, respectively. Compared to the model with ρ = 1 kN/mm 2 , the mid-span vertical deflections of the model with ρ = 10 kN/mm 2 on the 28th day and in the 3rd year are 13.08% smaller and 7.90% smaller, respectively; for the interface slip at the support, the vertical deflections are 89.13% smaller and 89.20% smaller, respectively; for the maximum compressive stress on the concrete slab, the vertical deflections are 1.97% larger and 1.91% larger, respectively; for the maximum tensile stress in the steel beam bottom flange, the vertical deflections are 0.96% smaller and 0.88% smaller, respectively.  = 0.5 kN/mm 2 on the 28th day and in the 3rd year are 13.90% larger and 8.40% larger, respectively; for the interface slip at the beam end, the vertical deflections are 91.64% larger and 90.62% larger, respectively; for the maximum compressive stress on the concrete slab, the vertical deflections are 2.2% smaller and 2.13% smaller, respectively; for the maximum tensile stress in the steel beam bottom flange, the vertical deflections are 1.06% larger and 0.96% larger, respectively. Compared to the model with  = 1 kN/mm 2 , the midspan vertical deflections of the model with  = 10 kN/mm 2 on the 28th day and in the 3rd year are 13.08% smaller and 7.90% smaller, respectively; for the interface slip at the support, the vertical deflections are 89.13% smaller and 89.20% smaller, respectively; for the maximum compressive stress on the concrete slab, the vertical deflections are 1.97% larger and 1.91% larger, respectively; for the maximum tensile stress in the steel beam bottom flange, the vertical deflections are 0.96% smaller and 0.88% smaller, respectively.

Time-Dependent Analysis
Time-dependent analysis of the above composite box beams focuses on the mechanical responses on the 28th day and in the 3rd year. Figures 13-18 exhibit the time-dependent responses of the composite beams, including the vertical deflections, interface slip, warping intensity function on the concrete slab due to shear deformation, warping intensity function at the bottom flange of the steel beam due to shear deformation, and stresses on the concrete slab and at the bottom flange of steel beam, respectively. Figures 13a-18a show the distribution of these responses on the 28th day and in the 3rd year, while Figures  13b-18b show the changes over time. Figure 13 shows the time-dependent variation in the vertical deflection of the composite beam. Figure 13a indicates that the long-term vertical deflection in the 3rd year is 47.01% larger than that on the 28th day. The time-dependent variation in the deflection at mid-span is shown in Figure 13b. It can be noted that the vertical deflection increases gradually over time, but in a narrower range.

Time-Dependent Analysis
Time-dependent analysis of the above composite box beams focuses on the mechanical responses on the 28th day and in the 3rd year. Figures 13-18 exhibit the timedependent responses of the composite beams, including the vertical deflections, interface slip, warping intensity function on the concrete slab due to shear deformation, warping intensity function at the bottom flange of the steel beam due to shear deformation, and stresses on the concrete slab and at the bottom flange of steel beam, respectively. Figures 13a, 14a, 15a, 16a, 17a and 18a show the distribution of these responses on the 28th day and in the 3rd year, while Figures 13b, 14b, 15b, 16b, 17b and 18b show the changes over time.

Time-Dependent Analysis
Time-dependent analysis of the above composite box beams focuses on the mechanical responses on the 28th day and in the 3rd year. Figures 13-18 exhibit the time-dependent responses of the composite beams, including the vertical deflections, interface slip, warping intensity function on the concrete slab due to shear deformation, warping intensity function at the bottom flange of the steel beam due to shear deformation, and stresses on the concrete slab and at the bottom flange of steel beam, respectively. Figures 13a-18a show the distribution of these responses on the 28th day and in the 3rd year, while Figures  13b-18b show the changes over time. Figure 13 shows the time-dependent variation in the vertical deflection of the composite beam. Figure 13a indicates that the long-term vertical deflection in the 3rd year is 47.01% larger than that on the 28th day. The time-dependent variation in the deflection at mid-span is shown in Figure 13b. It can be noted that the vertical deflection increases gradually over time, but in a narrower range.   Figure 14 shows the time-dependent interface slip between the concrete slab and the steel beam. This response is consistent with the common observation that obvious slip often occurs at the end of composite beams. The interface slip at the end of the composite beam continually decreases from 1.665 mm on the 28th day to 1.482 mm in the 3rd year, as illustrated in Figure 14a. Figure 14b shows the variation in the interface slip at the beam end. The interface slip gradually decreases over time, but in a narrower range.  Figure 15 shows the time-dependent variation in the shear lag warping intensity function of the concrete slab. Figure 15a shows that this response is more significant at the support of the composite beam, while this response is 0 at mid-span. This response in the 3rd year was 111.64% higher than that on the 28th day. Figure 15b exhibits the time-dependent variations in the warping intensity function of the concrete slab at the support of the composite beam. The warping intensity function of the concrete slab increases gradually over time, but in a narrow range. The increase slowed after the 400th day.   Figure 15 shows the time-dependent variation in the shear lag warping intensity function of the concrete slab. Figure 15a shows that this response is more significant at the support of the composite beam, while this response is 0 at mid-span. This response in the 3rd year was 111.64% higher than that on the 28th day. Figure 15b exhibits the time-dependent variations in the warping intensity function of the concrete slab at the support of the composite beam. The warping intensity function of the concrete slab increases gradually over time, but in a narrow range. The increase slowed after the 400th day.   Figure 16 shows the time-dependent variations in the warping intensity function at the bottom flange of the steel beam due to shear deformation. As shown in Figure 16a, this response at the beam end in the 3rd year is 7.01% higher than that on the 28th day. Figure 16b shows the changes in the response at the composite beam end. The response at the support of the composite beam gradually increases over time, but in a narrower range. The increase in the warping intensity function of the concrete slab is more significant than that at the bottom flange of the steel beam.  Figure 17 shows the stress distribution on the concrete slab and the variation with time. According to Figure 17a, the compressive stress at the edge of the concrete slab decreases from 5.829 MPa on the 28th day to 5.428 MPa in the 3rd year; the compressive stress at the intersection between the steel web and the concrete slab decreased from 6.075 time. According to Figure 17a, the compressive stress at the edge of the concrete slab decreases from 5.829 MPa on the 28th day to 5.428 MPa in the 3rd year; the compressive stress at the intersection between the steel web and the concrete slab decreased from 6.075 MPa on the 28th day to 5.665 MPa in the 3rd year; the compressive stress of the centerline of the concrete slab decreases from 5.446 MPa on the 28th day to 5.058 MPa in the 3rd year. Figure 17b shows the variation in these stresses with time and reveals that the stress on the concrete slab slowly decreases after the 400th day.   Figure 18b shows the variation in these stresses with time and reveals that the stress of steel slowly increases after the 400th day. To better show the time-dependent influence, the composite beam on the 3650th day was analyzed. Figure 19 shows the stress distribution on the concrete slab on the 28th and 3650th days. Figure 20 shows the stress distribution at the steel bottom on the 28th and 3650th days. The stress on the concrete slab decreases over this period (28th to 3650th day), and the stress at the bottom flange of the steel beam increases. These results indicate that the creep of concrete produces an unloading effect on the concrete slab.   Figure 13 shows the time-dependent variation in the vertical deflection of the composite beam. Figure 13a indicates that the long-term vertical deflection in the 3rd year is 47.01% larger than that on the 28th day. The time-dependent variation in the deflection at mid-span is shown in Figure 13b. It can be noted that the vertical deflection increases gradually over time, but in a narrower range. Figure 14 shows the time-dependent interface slip between the concrete slab and the steel beam. This response is consistent with the common observation that obvious slip often occurs at the end of composite beams. The interface slip at the end of the composite beam continually decreases from 1.665 mm on the 28th day to 1.482 mm in the 3rd year, as illustrated in Figure 14a. Figure 14b shows the variation in the interface slip at the beam end. The interface slip gradually decreases over time, but in a narrower range. Figure 15 shows the time-dependent variation in the shear lag warping intensity function of the concrete slab. Figure 15a shows that this response is more significant at the support of the composite beam, while this response is 0 at mid-span. This response in the 3rd year was 111.64% higher than that on the 28th day. Figure 15b exhibits the timedependent variations in the warping intensity function of the concrete slab at the support of the composite beam. The warping intensity function of the concrete slab increases gradually over time, but in a narrow range. The increase slowed after the 400th day. Figure 16 shows the time-dependent variations in the warping intensity function at the bottom flange of the steel beam due to shear deformation. As shown in Figure 16a, this response at the beam end in the 3rd year is 7.01% higher than that on the 28th day. Figure 16b shows the changes in the response at the composite beam end. The response at the support of the composite beam gradually increases over time, but in a narrower range. The increase in the warping intensity function of the concrete slab is more significant than that at the bottom flange of the steel beam. Figure 17 shows the stress distribution on the concrete slab and the variation with time. According to Figure 17a, the compressive stress at the edge of the concrete slab decreases from 5.829 MPa on the 28th day to 5.428 MPa in the 3rd year; the compressive stress at the intersection between the steel web and the concrete slab decreased from 6.075 MPa on the 28th day to 5.665 MPa in the 3rd year; the compressive stress of the centerline of the concrete slab decreases from 5.446 MPa on the 28th day to 5.058 MPa in the 3rd year. Figure 17b shows the variation in these stresses with time and reveals that the stress on the concrete slab slowly decreases after the 400th day. Figure 18 exhibits the stress distribution of the bottom flange of the steel beam and the change over time. Figure 18a shows that the stress at the steel web-concrete slab intersection increases from 96.349 MPa on the 28th day to 100.739 MPa in the 3rd year; the stress of the centerline of the steel bottom increases from 87.039 MPa on the 28th day to 91.133 MPa in the 3rd year. Figure 18b shows the variation in these stresses with time and reveals that the stress of steel slowly increases after the 400th day.
To better show the time-dependent influence, the composite beam on the 3650th day was analyzed. Figure 19 shows the stress distribution on the concrete slab on the 28th and 3650th days. Figure 20 shows the stress distribution at the steel bottom on the 28th and 3650th days. The stress on the concrete slab decreases over this period (28th to 3650th day), and the stress at the bottom flange of the steel beam increases. These results indicate that the creep of concrete produces an unloading effect on the concrete slab.  Figure 18a shows that the stress at the steel web-concrete slab intersection increases from 96.349 MPa on the 28th day to 100.739 MPa in the 3rd year; the stress of the centerline of the steel bottom increases from 87.039 MPa on the 28th day to 91.133 MPa in the 3rd year. Figure 18b shows the variation in these stresses with time and reveals that the stress of steel slowly increases after the 400th day. To better show the time-dependent influence, the composite beam on the 3650th day was analyzed. Figure 19 shows the stress distribution on the concrete slab on the 28th and 3650th days. Figure 20 shows the stress distribution at the steel bottom on the 28th and 3650th days. The stress on the concrete slab decreases over this period (28th to 3650th day), and the stress at the bottom flange of the steel beam increases. These results indicate that the creep of concrete produces an unloading effect on the concrete slab.

Conclusions
A beam finite element model with 18 DOFs considering slip, shear lag, and timedependent effects of the composite box beam is proposed in this study. The conclusions from the numerical results are drawn as follows: (1) The analytical model of the composite box beam is solved in the time domain by an

Conclusions
A beam finite element model with 18 DOFs considering slip, shear lag, and timedependent effects of the composite box beam is proposed in this study. The conclusions from the numerical results are drawn as follows: (1) The analytical model of the composite box beam is solved in the time domain by an accurate step-by-step method, and the Dirichlet series is employed into the creep function without the storage of stress and strain histories. In the spatial domain, the FEM discretizes the composite beam into finite two-node beam elements with 18 DOFs. An effective recursion method was developed to solve these equations of the system. (2) A fine numerical simulation by ANSYS is performed to validate the correctness and applicability of the proposed model in instantaneous analysis. Through a comparison between the classical long-term test results of the composite beams and the calculation results of the proposed model involving beam finite elements, the validity and applicability of the proposed model used in the long-term analysis was proven. (3) The proposed model is then applied to the analysis of the time-dependent behavior of composite beams. Some mechanical responses, including vertical deflections, interface slip, warping intensity function due to shear deformation, the stress on the concrete slab, and the stress of the steel beam, are analyzed. In this study, the variation in these responses between the 28th day and the 3rd year was the focus. It was determined that the shear connection stiffness, shrinkage, and creep all have a significant impact on these responses. (4) Compared with the model with ρ = 1 kN/mm 2 , the vertical deflection at mid-span in the 3rd year of the model with ρ = 0.5 kN/mm 2 is 8.40% larger; the interface slip at the end is 90.62% larger; the mid-span vertical deflection in the 3rd year of the model with ρ = 10 kN/mm 2 is 7.90% less; and the interface slip at the end is 89.20% smaller. (5) For the model with ρ = 1 kN/mm 2 , from the 28th day to the 3rd year, the mid-span vertical deflection increased by 47.01%, the interface slip at the end decreased by 10.99%, the warping intensity function of the concrete slab at end of the steel beam due to shear deformation increased by 111.64%, the warping intensity function at the steel bottom near the beam end increased by 7.01%, the maximum compressive stress on the concrete slab decreased by 6.75%, and the maximum tensile stress at the steel bottom flange increased by 4.56%. The above mechanical responses have been developing, and even in the 10th year, the development rate is still not negligible.