**Figure 5.** Buckling modes for the 300 mm-long C-section, SAM results. Subsequent modes are as follows: Local *m* > 1 (Mode1), primary global (Mode2), higher global distortional-flexural (Mode3), and third global distortional-flexural (Mode4). **Figure 5.** Buckling modes for the 300 mm-long C-section, SAM results. Subsequent modes are as follows: Local *m* > 1 (Mode1), primary global (Mode2), higher global distortional-flexural (Mode3), and third global distortional-flexural (Mode4). **Figure 5.** Buckling modes for the 300 mm-long C-section, SAM results. Subsequent modes are as follows: Local *m* > 1 (Mode1), primary global (Mode2), higher global distortional-flexural (Mode3), and third global distortional-flexural (Mode4).

**Figure 6.** Buckling modes for the 500 mm-long C-section, SAM results. The same notation of modes **Figure 6.** Buckling modes for the 500 mm-long C-section, SAM results. The same notation of modes as in Figure 5. **Figure 6.** Buckling modes for the 500 mm-long C-section, SAM results. The same notation of modes as in Figure 5.

The next stage consisted of an analysis of post-buckling equilibrium paths for the C-columns. Two-mode (*J* = 2 in Equation (5)) and three-mode (*J* = 3) approaches were considered for the SAM to facilitate an interpretation of the results. Moreover, the initial imperfections ζr \* = |0.2|, where the signs were selected in the most disadvantageous way [30–32,34], were assumed.

Figure 7 presents a collection of the SAM and FEM solutions showing post-buckling equilibrium paths in the system of the dimensionless load *P*/*P*min = *P*/*P*1 as a function of the dimensionless shortening ∆/∆*min* = ∆/∆1 (where *P*min corresponds to the lowest value of the bifurcation load, and ∆*min* denotes the minimal shortening caused by the load *P*min) for the 300 mm-long column. The post-buckling paths for a coupled interaction of two modes (*J* = 2): The local mode (*i* = 1) with the primary global one (*i* = 2) (denoted here as SAM-1,2) and for two cases of interaction of three modes (*J* = 3): The local mode (*i* = 1) with the two lowest global modes (i.e., for *i* = 2,3) (denoted as SAM-1,2,3), and the local mode (*i* = 1) with the global primary mode (*i* = 2) and the third mode (*i* = 4) (denoted as SAM-1,2,4), are presented. For the two-mode approach, the post-buckling path grows monotonically. For the three-mode approach, the SAM-1,2,3 curve flattens, not reaching the ultimate load-carrying capacity, whereas, for the SAM-1,2,4, it attains the ultimate value of *P*/*P*min ≈ 1.1 and then falls, which corresponds

C-section.

C-section.

to an unstable equilibrium path. It follows from these comparisons that an effect of the third global mode (*i* = 4) is crucial for interactive buckling of the 300 mm-long C-column. An application of the SAM thus enables an easy analysis of the effect of the modes under study on post-buckling paths, which is not possible virtually with the FEM. The FEM was used to verify the SAM results. In the FEM calculations, the amplitude of the initial deflection equal to 0.2 of the component plate thickness was assumed for the local mode (*i* = 1) with the global primary mode (*i* = 2) and the second mode (*i* = 3) or third global mode (*i* = 4). The results were the same. A very good agreement between the FEM curve and the SAM-1,2,3 curve was obtained. The post-buckling curve SAM-1,2 begins to differ more and more from the FEM curve for the load *P*/*P*min ≥ 1.8. It was not possible to verify the curve SAM-1,2,4 with the FEM. Attention should be paid to the FEM curve. It breaks catastrophically for *P*/*P*min ≈ 1.95 and ∆/∆*min*≈ 3.4 due to a lack of convergence of the Riks method. Thus, it is the ultimate value. None of the changes in the parameters of the Riks method allowed the calculations to continue. The authors put forward a hypothesis that when the secondary (*i* = 3) and third (*i* = 4) global modes begin to interact, numerical convergence of the FEM is lost, and the computations are catastrophically disrupted. This hypothesis requires further thorough investigations to be carried out. In Figure 8, a mode for the ultimate FEM value is shown. As can be easily seen in the central part of the C-channel, deflections of the plates and the web flattened out, with the symmetry of displacements remaining with respect to the cross-section axis. *Materials* **2020**, *13*, x 10 of 17

**Figure 7.** Post-buckling equilibrium paths in the dimensionless system for the 300 mm-long C-section.

**Figure 8.** Mode of the 300-mm long C-section in the ultimate state, FEM results. **Figure 8.** Mode of the 300-mm long C-section in the ultimate state, FEM results.

In the next figures (Figures 9 and 10), results for the 500 mm-long C-channel are presented. In Figure 9, the post-buckling equilibrium paths *P*/*P*min = *P*/*P*1 as a function of the shortening ∆/∆*min* = ∆/∆1 are shown for the SAM and the FEM. The curves SAM-1,2 for the two-mode approach (*J* = 2) and SAM-1,2,3 for the three-mode approach (*J* = 3) are situated close to each other up to the load *P*/*P*min ≈ 2.5. The post-buckling curve SAM-1,2,4 differs from the remaining curves already at *P*/*P*min > 0.85 and attains the ultimate load-carrying capacity for the load *P*/*P*min ≈ 2 and the shortening ∆/∆*min* ≈ 5.1. The FEM post-buckling equilibrium path is similar to the curves SAM-1,2 and SAM-1,2,3, but lies above them. In this case, we were not able to obtain a curve close to SAM-1,2,4 either. For example, for the 300 mm-long C-column, the FEM computations were disrupted catastrophically by a lack of numerical convergence for *P*/*P*min ≈ 1.35 and ∆/∆*min*≈ 1.65, which should be treated as the ultimate value. Figure 10 presents the buckling mode of the C-section at the ultimate point.

1

**Figure 10.** Mode of the 500 mm-long C-section in the ultimate state, FEM results.

C-section.

Non-symmetrical buckling of C-section strips, which can follow from an interaction of anti-symmetrical modes, can be observed. They have been neglected in the present study.

**Figure 8.** Mode of the 300-mm long C-section in the ultimate state, FEM results.

**Figure 7.** Post-buckling equilibrium paths in the dimensionless system for the 300 mm-long

*Materials* **2020**, *13*, x 10 of 17

**Figure 9.** Post-buckling equilibrium paths in the dimensionless system for the 500 mm-long C-section. C-section.

**Figure 10.** Mode of the 500 mm-long C-section in the ultimate state, FEM results. **Figure 10.** Mode of the 500 mm-long C-section in the ultimate state, FEM results.

For both the lengths of C-columns, a huge impact of high values of global buckling modes, which are unrealistic in the theory of interactive buckling from the practical point of view, can be seen. The present study is just an attempt to identify and to describe the phenomenon of the influence of secondary and third global distortional-flexural modes on the load-carrying capacity of medium-long C-channels.

#### *3.2. TH-Channel Columns*

1 In this section, the stability and post-buckling equilibrium paths of "top-hat"- section columns (Figure 1b) of the following lengths: 300 and 500 mm, will be discussed. Table 1 includes the bifurcation values of stresses obtained with the SAM and the FEM. Moreover, for the local mode (*i* = 1), the number of halfwaves m along the longitudinal direction is given in brackets. The same index notations are used as for the C-channels.

A very good agreement of the SAM and FEM results was obtained, similarly to the case of the C-channels. It should be noticed that an introduction of edge reinforcements in TH-section columns in comparison to C-section columns results is almost twice as high an increase in local bifurcational loads. What is more, the values of primary global loads for both lengths are very close. The ratio of the value of the second global mode to the local one is more than seven times higher, whereas for the third mode, it is 40 times higher, correspondingly.

Figure 11 shows bifurcation modes for the 300 mm column, whereas Figure 12 shows those for 500 mm. Local modes (Figures 11a and 12a) are the 'pure' local ones, for which the edges do not displace. They differ in the number of halfwaves along the longitudinal direction. The amplitudes of modes (i.e., the shape of the mode in the cross-section) are identical and, analogously, primary global modes (*i* = 2) (Figures 11b and 12b). In this case, the mode is distortional as the reinforced

edge corner rotates together with the plate with respect to the web. Global secondary modes (*i* = 3) (Figures 11c and 12c) are the distortional-flexural ones and they are similar, provided the interpretation for the C-channel from Figure 3a,b is assumed. The global secondary mode (*i* = 3) for the 500 mm column is even "purer" than that for 300 mm, i.e., it is closer to the flexural one. Third modes (*i* = 4) are the distortional-flexural ones. In Figures 13 and 14, buckling modes of columns 300 and 500 mm long, respectively, are presented for the SAM. When Figures 11–14 are compared, it can be easily noticed that the buckling modes determined with the SAM and FEM agree very well. *Materials* **2020**, *13*, x 12 of 17 *Materials* **2020**, *13*, x 12 of 17

**Figure 11.** Buckling modes for the 300 mm-long TH-channel column, FEM results. Subsequent modes are as follows: Local *m* > 1 (**a**), primary global (**b**), higher global distortional-flexural (**c**), and third global distortional-flexural (**d**). **Figure 11.** Buckling modes for the 300 mm-long TH-channel column, FEM results. Subsequent modes are as follows: Local *m* > 1 (**a**), primary global (**b**), higher global distortional-flexural (**c**), and third global distortional-flexural (**d**). modes are as follows: Local *m* > 1 (**a**), primary global (**b**), higher global distortional-flexural (**c**), and third global distortional-flexural (**d**).

(**c**) (**d**) **Figure 12.** Buckling modes for the 500 mm-long TH-channel column, FEM results. Subsequent modes are as follows: Local *m* > 1 (**a**), primary global (**b**), higher global distortional-flexural (**c**), and third global distortional-flexural (**d**). **Figure 12.** Buckling modes for the 500 mm-long TH-channel column, FEM results. Subsequent modes are as follows: Local *m* > 1 (**a**), primary global (**b**), higher global distortional-flexural (**c**), and third global distortional-flexural (**d**). **Figure 12.** Buckling modes for the 500 mm-long TH-channel column, FEM results. Subsequent modes are as follows: Local *m* > 1 (**a**), primary global (**b**), higher global distortional-flexural (**c**), and third global distortional-flexural (**d**).

**Figure 13.** Buckling modes for the 300 mm-long TH-channel column, SAM results. The same notation of modes as in Figure 5. **Figure 13.** Buckling modes for the 300 mm-long TH-channel column, SAM results. The same notation of modes as in Figure 5. **Figure 13.** Buckling modes for the 300 mm-long TH-channel column, SAM results. The same notation of modes as in Figure 5.

**Figure 14.** Buckling modes for the 500 mm-long TH-channel column, SAM results. The same notation of modes as in Figure 5. **Figure 14.** Buckling modes for the 500 mm-long TH-channel column, SAM results. The same notation of modes as in Figure 5. **Figure 14.** Buckling modes for the 500 mm-long TH-channel column, SAM results. The same notation of modes as in Figure 5.

The SAM and FEM calculations for the post-buckling states were conducted analogously as for the C-channels, that is to say, the same level of mode imperfections and two- and three-mode approaches for the SAM were assumed. The post-buckling paths in the systems *P/P*min = *P/P*1 versus the shortening *Δ/Δmin = Δ/Δ1* are presented in Figure 15 for 300 mm and in Figure 16 for 500 mm. For both the lengths, the ultimate load-carrying capacity was obtained both for the two-mode (*J =* 2) and three-mode (*J =* 3) approaches. The curves SAM-1,2, SAM-1,2,3, and SAM-1,2,4 overlap up to the ultimate point, and then they fall, but at 300 mm, the SAM-1,2,4 (Figure 15) falls most dramatically. For the length of 500 mm, the curves SAM-1,2,3 and SAM-1,2,4 fall much more sharply than SAM-1,2. For 300 mm (Figure 15), the ultimate point occurs at *P/P*min ≈ 1.7 and Δ/Δmin ≈ 2.25, whereas for 500 mm (Figure 16), it occurs at *P/P*min ≈ 2 and Δ/Δmin ≈ 2.6. The FEM paths obtained for 300 (Figure 15) and 500 mm (Figure 16) are close to the ultimate load-carrying capacity for the SAM, which occurs at *P/P*min > 2.5 and *Δ/Δmin* > 5.4. Similarly, as for the C-channels, the FEM computations were catastrophically disrupted due to a lack of numerical convergence. The SAM and FEM calculations for the post-buckling states were conducted analogously as for the C-channels, that is to say, the same level of mode imperfections and two- and three-mode approaches for the SAM were assumed. The post-buckling paths in the systems *P/P*min = *P/P*1 versus the shortening *Δ/Δmin = Δ/Δ1* are presented in Figure 15 for 300 mm and in Figure 16 for 500 mm. For both the lengths, the ultimate load-carrying capacity was obtained both for the two-mode (*J =* 2) and three-mode (*J =* 3) approaches. The curves SAM-1,2, SAM-1,2,3, and SAM-1,2,4 overlap up to the ultimate point, and then they fall, but at 300 mm, the SAM-1,2,4 (Figure 15) falls most dramatically. For the length of 500 mm, the curves SAM-1,2,3 and SAM-1,2,4 fall much more sharply than SAM-1,2. For 300 mm (Figure 15), the ultimate point occurs at *P/P*min ≈ 1.7 and Δ/Δmin ≈ 2.25, whereas for 500 mm (Figure 16), it occurs at *P/P*min ≈ 2 and Δ/Δmin ≈ 2.6. The FEM paths obtained for 300 (Figure 15) and 500 mm (Figure 16) are close to the ultimate load-carrying capacity for the SAM, which occurs at *P/P*min > 2.5 and *Δ/Δmin* > 5.4. Similarly, as for the C-channels, the FEM computations were catastrophically disrupted due to a lack of numerical convergence. The SAM and FEM calculations for the post-buckling states were conducted analogously as for the C-channels, that is to say, the same level of mode imperfections and two- and three-mode approaches for the SAM were assumed. The post-buckling paths in the systems *P*/*P*min = *P*/*P*1 versus the shortening ∆/∆*min* = ∆/∆1 are presented in Figure 15 for 300 mm and in Figure 16 for 500 mm. For both the lengths, the ultimate load-carrying capacity was obtained both for the two-mode (*J* = 2) and three-mode (*J* = 3) approaches. The curves SAM-1,2, SAM-1,2,3, and SAM-1,2,4 overlap up to the ultimate point, and then they fall, but at 300 mm, the SAM-1,2,4 (Figure 15) falls most dramatically. For the length of 500 mm, the curves SAM-1,2,3 and SAM-1,2,4 fall much more sharply than SAM-1,2. For 300 mm (Figure 15), the ultimate point occurs at *P*/*P*min ≈ 1.7 and ∆/∆min ≈ 2.25, whereas for 500 mm (Figure 16), it occurs at *P*/*P*min ≈ 2 and ∆/∆min ≈ 2.6. The FEM paths obtained for 300 (Figure 15) and 500 mm (Figure 16) are close to the ultimate load-carrying capacity for the SAM, which occurs at *P*/*P*min > 2.5 and ∆/∆*min* > 5.4. Similarly, as for the C-channels, the FEM computations were catastrophically disrupted due to a lack of numerical convergence.

**Figure 15.** Post-buckling equilibrium paths in the dimensionless system for the 300 mm-long TH-section.

TH-section.

TH-section.

**Figure 15.** Post-buckling equilibrium paths in the dimensionless system for the 300 mm-long

*Materials* **2020**, *13*, x 14 of 17

**Figure 16.** Post-buckling equilibrium paths in the dimensionless system for the 500 mm-long TH-section. mm C-channel (Figure 10). On the other hand, the mode for 500 mm (Figure 18) is symmetrical in practice.

In Figures 17 and 18, buckling modes for the FEM ultimate points at 300 and 500 mm, respectively, are presented. The mode for 300 mm (Figure 17) is more asymmetrical than for the 500 mm C-channel (Figure 10). On the other hand, the mode for 500 mm (Figure 18) is symmetrical in practice. In the FEM, the modes considered in the analysis are not under control as opposed to the SAM. The asymmetrical mode (Figure 17) could be attained via an automatic introduction of anti-symmetrical eigenvalues (i.e., bifurcation ones) into the computations. Edge reinforcements for the top-hat channels result in higher load-carrying capacity when compared to the C-channels.

**Figure 17.** Mode of the 300 mm-long TH-channel column in the ultimate state, FEM results. **Figure 17.** Mode of the 300 mm-long TH-channel column in the ultimate state, FEM results. *Materials* **2020**, *13*, x 15 of 17

**Figure 18.** Mode of the 500 mm-long TH-channel column in the ultimate state, FEM results. **Figure 18.** Mode of the 500 mm-long TH-channel column in the ultimate state, FEM results.

In the SAM calculations for C- and TH-channels for the SAM-1,2,4 three-mode approach, an unexpected effect of the third global distortional-flexural mode on the interactive buckling and the load-carrying capacity was found. It was not possible to verify the SAM calculations within the FEM; nevertheless, they unexpectedly cast light on a catastrophic disruption of the computations In the FEM, the modes considered in the analysis are not under control as opposed to the SAM. The asymmetrical mode (Figure 17) could be attained via an automatic introduction of anti-symmetrical eigenvalues (i.e., bifurcation ones) into the computations. Edge reinforcements for the top-hat channels result in higher load-carrying capacity when compared to the C-channels.

for post-buckling paths due to a lack of convergence, although Rik's procedure with various initial parameters was applied. The observed numerical phenomena for the SAM and the FEM thus require further intensive investigations with respect to modeling and the methods applied. **4. Conclusions**  Multimodal buckling of the C- and TH-section (the so-called "top-hat") columns, made of isotropic material and subjected to compression, was considered. The column lengths were selected In the SAM calculations for C- and TH-channels for the SAM-1,2,4 three-mode approach, an unexpected effect of the third global distortional-flexural mode on the interactive buckling and the load-carrying capacity was found. It was not possible to verify the SAM calculations within the FEM; nevertheless, they unexpectedly cast light on a catastrophic disruption of the computations for post-buckling paths due to a lack of convergence, although Rik's procedure with various initial parameters was applied. The observed numerical phenomena for the SAM and the FEM thus require further intensive investigations with respect to modeling and the methods applied.

equilibrium paths prepared within the SAM in the case of an interaction of the lowest local and other modes, when higher global distortional-flexural modes are neglected, the curves grow monotonically. However, if the global distortional flexural buckling mode is considered, the post-buckling equilibrium path grows up to the ultimate point, and then, it falls. This regularity repeats in all the columns under analysis. Moreover, in both the methods, attention is drawn to their conditions, particularly to the FEM catastrophic numerical behavior in all the cases of computations under discussion. Hence, the observed phenomena require further thorough and complex investigations, including the experimental ones. Finally, it should be added that all the considerations presented here refer to the elastic state. Thus, all plastic effects and plastic

**Author Contributions:** Conceptualization, Z.K.; methodology, A.T. and Z.K.; software, A.T. and Z.K.; formal analysis, A.T. and Z.K.; writing—original draft, A.T. and Z.K.; writing—review and editing, A.T. All authors

**Funding:** The studies have been conducted as a part of the research projects financed by the National Science

1. Van der Heijden, A.M.A. *Koiter's Elastic Stability of Solids and Structures*; Publisher: Cambridge Universitet

2. Doyle, J.F. *Nonlinear Analysis of Thin-Walled Structures*; Springer-Verlag: New York, USA, 2001,

destruction mechanisms that occurred have been neglected.

have read and agree to the published version of the manuscript.

Centre, Poland - decision numbers: UMO-2017/25/B/ST8/01046. **Conflicts of Interest:** The authors declare no conflict of interest

doi:10.1007/978-1-4757-3546-8.

Delft, Delft, The Netherlands, 2008, doi:10.1017/CBO9780511546174.

**References**

to observe strong interactions between various eigenvalues. The problem of nonlinear stability was solved with two methods, namely the semi-analytical method (SAM) using Koiter's perturbation

## **4. Conclusions**

Multimodal buckling of the C- and TH-section (the so-called "top-hat") columns, made of isotropic material and subjected to compression, was considered. The column lengths were selected to observe strong interactions between various eigenvalues. The problem of nonlinear stability was solved with two methods, namely the semi-analytical method (SAM) using Koiter's perturbation approach and the finite element method (FEM). The results discussed prove that there is a strong effect of higher global distortional-flexural buckling modes (i.e., secondary and third) on the interactive buckling of Cand TH-columns and their load-carrying capacity. On the post-buckling equilibrium paths prepared within the SAM in the case of an interaction of the lowest local and other modes, when higher global distortional-flexural modes are neglected, the curves grow monotonically. However, if the global distortional flexural buckling mode is considered, the post-buckling equilibrium path grows up to the ultimate point, and then, it falls. This regularity repeats in all the columns under analysis. Moreover, in both the methods, attention is drawn to their conditions, particularly to the FEM catastrophic numerical behavior in all the cases of computations under discussion. Hence, the observed phenomena require further thorough and complex investigations, including the experimental ones. Finally, it should be added that all the considerations presented here refer to the elastic state. Thus, all plastic effects and plastic destruction mechanisms that occurred have been neglected.

**Author Contributions:** Conceptualization, Z.K.; methodology, A.T. and Z.K.; software, A.T. and Z.K.; formal analysis, A.T. and Z.K.; writing—original draft, A.T. and Z.K.; writing—review and editing, A.T. All authors have read and agree to the published version of the manuscript.

**Funding:** The studies have been conducted as a part of the research projects financed by the National Science Centre, Poland—decision numbers: UMO-2017/25/B/ST8/01046.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Finite Beam Element for Curved Steel–Concrete Composite Box Beams Considering Time-Dependent E**ff**ect**

## **Guang-Ming Wang, Li Zhu \* , Xin-Lin Ji and Wen-Yu Ji**

School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China<EMAIL_ADDRESS>(G.-M.W<EMAIL_ADDRESS>(X.-L.J<EMAIL_ADDRESS>(W.-Y.J.)

**\*** Correspondence<EMAIL_ADDRESS>
Received: 19 June 2020; Accepted: 20 July 2020; Published: 22 July 2020

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

**Keywords:** curved steel–concrete composite box beam; two-node finite beam element with 26 DOFs; long-term behavior; age-adjusted effective modulus method

## **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–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. 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.

straight composite beams. Due to the complex mechanical behavior and time-dependent effect of

*Materials* **2020**, *12*, x FOR PEER REVIEW 3 of 17

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

**Figure 1.** Geometric dimensions and coordinate system of the curved composite box beam: (**a**) cross section, (**b**) plane. **Figure 1.** Geometric dimensions and coordinate system of the curved composite box beam: (**a**) cross section, (**b**) plane.

#### *2.3. Displacement Modes and Strain Components of Curved Composite Box Beams 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 *ak*Ω*x*(*z*) and *ak*Ω*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, *fk* (*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 *Wzk*(*z*,*s*), tangential displacement *Wsk*(*z*,*s*), and normal displacement *Wnk*(*z*,*s*) at any point of the composite box beam can be obtained as follows:

$$\begin{aligned} \mathcal{W}\_{zk}(\mathbf{z}, \mathbf{s}) &= w - (u' + wk\_0 + a\_k \Omega\_z k\_0)\mathbf{x} - v'y + a\_k \Omega\_z \\ -\omega(\mathbf{x}, y)(\boldsymbol{\beta} + v'k\_0) - \omega\_\mathbf{d}(\mathbf{x}, y)\theta\_\mathbf{d}\mathbf{'} + \psi\_k(\mathbf{x})f\_k \end{aligned} \tag{1}$$

$$\mathcal{W}\_{\rm sk}(z,s) = -(u + a\_{\rm k}\Omega\_{\rm x})\sin\alpha - v\cos\alpha + \theta\rho\_{\rm s} + \theta\_{\rm d}D\_{\rm s} \tag{2}$$

$$W\_{nk}(z,s) = (u + a\_k \Omega\_{\mathbf{x}}) \cos \alpha - v \sin \alpha - \theta \rho\_n - \theta\_{\mathbf{d}} D\_n \tag{3}$$

where *k*0 = 1/*R*0 is the initial curvature of the centroid line of the box beam sections; *R*0 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), ρ*s* = (*y* − *y*s) sin α − (*x* − *x*s) cos α is the tangential displacement distribution of the cross-sections of the composite beam under a unit torsion angle; ρ*n* = (*y* − *y*s) cos α + (*x* − *x*s) sin α is the normal displacement distribution of the cross-sections of the composite beam under a unit torsion angle; *Ds* = [(*y* − *y*d) sin α − (*x* − *x*d) cos α]Ψd is the tangential displacement distribution of the composite beam under a unit distortion angle; *Dn* = [(*y* − *y*d) cos α + (*x* − *x*d) sin α]Ψ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; Ψd = −1/(1 + ν) is the distortion parameter of the top and bottom plates, and Ψd = ν/(1 + ν) is the distortion parameter of the web, in which ν = (2*h*1*b*t + 2*h*3*b*s)/(*h*2*c* + *h*4*c*) is a constant related only to the section shape [20]; *h*1 and *h*3 are the vertical distances from the distortion center to the top and bottom plates, respectively; and *h*2 and *h*4 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:

$$
\mu\_{\rm sp} = \mathcal{D}\Omega\_{\rm x} \tag{4}
$$

$$w\_{\rm sp} = \mathcal{D}\Omega\_{\rm z} \tag{5}$$

where *u*sp and *w*sp denote the transverse and longitudinal interfacial slip, respectively.

*Materials* **2020**, *13*, 3253

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}{lcl} \varepsilon\_{\rm pk} &= \frac{\partial W\_{\rm x}}{\rho(\mathbf{x})\partial\phi} + \frac{W\_{\rm s}\sin\alpha}{\rho(\mathbf{x})} - \frac{W\_{\rm n}\cos\alpha}{\rho(\mathbf{x})} \\ &= \frac{R\_{0}}{\rho(\mathbf{x})} \left( \frac{\partial W\_{\rm x}}{R\_{0}\partial\phi} + \frac{W\_{\rm s}\sin\alpha}{R\_{0}} - \frac{W\_{\rm n}\cos\alpha}{R\_{0}} \right) \\ &= \frac{R\_{0}}{\rho(\mathbf{x})} \left[ \mathbf{w}^{\prime} - \mathbf{x}\mathbf{u}^{\prime\prime} - \mathbf{y}\mathbf{v}^{\prime\prime} + a\_{\rm k}(1 - k\_{0}\mathbf{x})\Omega\_{\rm z}{}^{\prime} - k\_{0}u - k\_{0}\mathbf{x}\mathbf{w}^{\prime} - a\_{\rm k}k\_{0}\Omega\_{\rm x} + \\ & \qquad \qquad \qquad k\_{0}(\mathbf{y} - \mathbf{y}\_{\rm s})\Theta - \omega(\mathbf{y}^{\prime} + k\_{0}\mathbf{z}^{\prime\prime}) + k\_{0}(\mathbf{y} - \mathbf{y}\_{\rm d})\Psi\_{\rm d}\Theta\_{\rm d} - \omega\_{\rm d}\mathbf{d}\_{\rm d}^{\prime\prime} + \psi\_{\rm f}\mathbf{f}\_{\rm s}^{\prime} \right] \end{array} \tag{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\_{\rm pk} = \frac{R\_0}{\rho(x)} \left[ (\theta' + \upsilon' k\_0) r^\* - (\beta + \upsilon' k\_0) \frac{\partial \omega}{\partial s} + \psi\_{k,x} f\_k \right] \tag{7}$$

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

$$\frac{\partial \boldsymbol{\omega}}{\partial \mathbf{s}} = \boldsymbol{r}^\* - \frac{\boldsymbol{\Omega}}{\mathbf{f} \boldsymbol{\Phi} (\mathbf{ds}/\mathbf{f})} \tag{8}$$

where Ω = H ˜*t*d*s* is twice the surrounding area, and ˜*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:

$$
\delta\!\!\!\/ \Pi = \underset{\begin{subarray}{c}\Phi\text{ }\delta\text{s}\_{\text{s}}^{\text{T}}\sigma\_{\text{s}}\text{d}\mu\rho(\mathbf{x})\text{d}\phi + \int\limits\_{\Phi} \underset{\begin{subarray}{c}\Phi\text{ }\delta\_{\text{c}}}{}\sigma\_{\text{c}}\text{d}\mu\rho(\mathbf{x})\text{d}\phi + \int\limits\_{\Phi} \underset{\begin{subarray}{c}\Phi\text{ }\delta\text{s}\_{\text{s}}\end{subarray}}{}\delta\mathbf{d}\_{\text{s}\text{l}\mathbf{s}}^{\text{T}}\mathbf{q}\_{\text{s}\text{h}}\mathbf{dx}\rho(\mathbf{x})\text{d}\phi + \\ \int\limits\_{\Phi} \delta\theta\_{\text{d}}\text{K}\_{\text{R}}\theta\_{\text{d}}\rho(\mathbf{x})\text{d}\phi - \sum\delta\mathbf{W}^{\text{T}}\mathbf{Q} - \int\limits\_{\Phi} \delta\mathbf{W}^{\text{T}}\mathbf{q}\rho(\mathbf{x})\text{d}\phi = 0 \\ \forall \delta\epsilon\_{\text{s}},\delta\epsilon\_{\text{c}},\delta d\_{\text{s}\text{l}\mathbf{p}},\delta\theta\_{\text{d}},\delta\mathbf{W}
\end{subarray}} (9)$$

where R Φ v *A*s δε T **s** σ**s**d*a*ρ(*x*)dφ and R Φ v *A*c δε T **c** σ**c**d*a*ρ(*x*)dφ 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:

$$
\varepsilon\_k = \mathbf{S} \mathbf{B}\_k \mathbf{d} \tag{10}
$$

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

$$\mathbf{S} = \begin{bmatrix} \frac{R\_0}{\rho(\mathbf{x})} & \frac{R\_0}{\rho(\mathbf{x})} \mathbf{x} & \frac{R\_0}{\rho(\mathbf{x})} \mathbf{y} & \frac{R\_0}{\rho(\mathbf{x})} \omega & \frac{R\_0}{\rho(\mathbf{x})} \omega\_\mathbf{d} & \frac{R\_0}{\rho(\mathbf{x})} \psi\_\mathbf{c} & \frac{R\_0}{\rho(\mathbf{x})} \psi\_\mathbf{s} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{R\_0}{\rho(\mathbf{x})} r^\* & \frac{\partial \omega}{\partial \mathbf{s}} & \psi\_{\mathbf{c},\mathbf{x}} & \psi\_{\mathbf{s},\mathbf{x}} \end{bmatrix} \tag{11}$$

$$d = \begin{array}{ccccc} \begin{array}{c} [u] \ \end{array} & [v] \end{array} \begin{array}{c} [w] \end{array} \quad \begin{array}{c} [\theta] \end{array} \quad \begin{array}{c} [\theta] \end{array} \quad \begin{array}{c} [\theta\_d] \end{array} \quad \begin{array}{c} [\Omega] \end{array} \quad \begin{array}{c} [f] \end{array} \end{array} \tag{12}$$

where

$$\begin{bmatrix} \boldsymbol{u} \end{bmatrix} = \begin{pmatrix} \boldsymbol{u} & \boldsymbol{u}' & \boldsymbol{u}'' \end{pmatrix} \tag{13}$$

$$\begin{bmatrix} \mathbf{v} \end{bmatrix} = \begin{pmatrix} \ v & v' & v'' \end{pmatrix} \tag{14}$$

$$\begin{bmatrix} \mathbf{w} \end{bmatrix} = \begin{pmatrix} \ w & w' \ \end{pmatrix} \tag{15}$$

$$\begin{bmatrix} \theta \end{bmatrix} = \begin{pmatrix} \ \theta & \theta' \ \end{pmatrix} \tag{16}$$

$$\begin{bmatrix} \mathfrak{g} \end{bmatrix} = \left( \begin{array}{cc} \beta & \beta' \end{array} \right) \tag{17}$$

$$\begin{bmatrix} \theta\_{\mathbf{d}} \end{bmatrix} = \left( \begin{array}{cc} \theta\_{\mathbf{d}} & \theta\_{\mathbf{d}}' & \theta\_{\mathbf{d}}'' \end{array} \right) \tag{18}$$

$$\begin{bmatrix} \mathbf{[0]} \end{bmatrix} = \left( \begin{array}{cc} \mathfrak{Q}\mathbf{1}\_{\mathbf{x}} & \mathfrak{Q}\mathbf{1}\_{\mathbf{x}} \end{array} \begin{array}{cc} \mathfrak{Q}\mathbf{1}\_{\mathbf{s}} \end{array} \begin{array}{cc} \mathfrak{Q}\mathbf{1}\_{\mathbf{s}} & \mathfrak{Q}\mathbf{1}\_{\mathbf{s}} \end{array} \right) \tag{19}$$

$$[\mathbf{f}] = \begin{pmatrix} f\_{\mathbf{c}} & f\_{\mathbf{c}}' & f\_{\mathbf{s}} & f\_{\mathbf{s}}' \end{pmatrix} \tag{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:

$$
\sigma\_{\rm s} = \left\{ \begin{array}{c} \sigma\_{\rm s} \\ \tau\_{\rm s} \end{array} \right\} = E\_{\rm s} \left[ \begin{array}{cc} 1 & 0 \\ 0 & 1/2(1+\mu\_{\rm s}) \end{array} \right] \left\{ \begin{array}{c} \varepsilon\_{\rm s} \\ \gamma\_{\rm s} \end{array} \right\} \tag{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} \varepsilon\_{\varepsilon}(t) - \varepsilon\_{\mathrm{c}}^{\mathrm{sh}}(t) \\ \gamma\_{\varepsilon}(t) \end{array} \right\} = J(t, t\_{0}) \begin{bmatrix} 1 & 0 \\ 0 & 2(1 + \mu\_{\mathrm{c}}) \end{bmatrix} \left\{ \begin{array}{c} \sigma\_{\varepsilon}(t\_{0}) \\ \tau\_{\mathrm{c}}(t\_{0}) \end{array} \right\} + \int\_{t\_{0}}^{t} J(t, \tau) \begin{bmatrix} 1 & 0 \\ 0 & 2(1 + \mu\_{\mathrm{c}}) \end{bmatrix} \left\{ \begin{array}{c} \mathrm{d}\sigma\_{\varepsilon}(\tau) \\ \mathrm{d}\tau\_{\mathrm{c}}(\tau) \end{array} \right\} \tag{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:

$$\sigma\_{\mathbf{c}}(t) = \begin{Bmatrix} \sigma\_{\mathbf{c}}(t) \\ \tau\_{\mathbf{c}}(t) \end{Bmatrix} = E\_{\mathbf{c}\mathbf{c}}(t, t\_0) \begin{bmatrix} 1 & 0 \\ 0 & 1/(2 + 2\mu\_{\mathbf{c}}) \end{bmatrix} \begin{Bmatrix} \varepsilon\_{\mathbf{c}}(t) - \varepsilon\_{\mathbf{c}}^{\mathrm{sh}}(t) \\ \gamma\_{\mathbf{c}}(t) \end{Bmatrix} + \xi(t, t\_0) \begin{Bmatrix} \sigma\_{\mathbf{c}}(t\_0) \\ \tau\_{\mathbf{c}}(t\_0) \end{Bmatrix} \tag{23}$$

where *t*0 is the initial loading age; ε sh c is the concrete shrinkage strain; µc is the Poisson's ratio of concrete; *J*(*t*,*t*0) 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*0) is the concrete age-adjusted effective elastic modulus; and ξ(*t*,*t*0) is the influence factor of the creep effects. *E*ce(*t*,*t*0) and ξ(*t*,*t*0) can be calculated by Equations (24) and (25), respectively.

$$E\_{\rm ce}(t\_{\prime}t\_{0}) = \frac{E\_{\rm c}(t\_{0})}{1 + \chi(t\_{\prime}t\_{0})\varphi(t\_{\prime}t\_{0})} \tag{24}$$

$$
\xi(t, t\_0) = \frac{\varphi(t, t\_0)(\chi(t, t\_0) - 1)}{1 + \chi(t, t\_0)\varphi(t, t\_0)} \tag{25}
$$

where *E*c(*t*0) is the concrete elastic modulus at time *t*0, ϕ(*t*,*t*0) is the concrete creep coefficient, and χ(*t*,*t*0) is the concrete aging coefficient, which can be derived from Equation (26).

$$\chi(t, t\_0) = \frac{E\_\mathbf{c}(t\_0)}{E\_\mathbf{c}(t\_0) - \mathbf{R}(t, t\_0)} - \frac{1}{\varphi(t, t\_0)}\tag{26}$$

where *R*(*t*,*t*0) is the concrete relaxation function, which can be calculated according to the literature [21]. Note that R Φ R 2*b*ts δ*d* T **slip***q***sh**d*x*ρ(*x*)dφ 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:

$$d\_{\rm slip} = \left\{ \begin{array}{cc} u\_{\rm sp} & w\_{\rm sp} \end{array} \right\}^{\rm T} = \mathbf{S}\_{\rm \Omega} \mathbf{B}\_{\rm \Omega} d \tag{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:

$$\boldsymbol{q}\_{\mathbf{sh}} = \begin{Bmatrix} \boldsymbol{q}\_{\text{us}} & \boldsymbol{q}\_{\text{ws}} \end{Bmatrix}^{\mathrm{T}} = \begin{bmatrix} \rho\_{\text{u}} & 0 \\ 0 & \rho\_{\text{w}} \end{bmatrix} \begin{Bmatrix} \boldsymbol{u}\_{\text{sp}} \\ \boldsymbol{w}\_{\text{sp}} \end{Bmatrix} = \rho\_{\text{sllip}} d\_{\text{sllip}} \tag{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 R Φ δθd*K*Rθdρ(*x*)dφ 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 P δ*W*T*Q* and R Φ δ*W*T*q*ρ(*x*)dφ in Equation (9) are the external virtual work produced

by external loads, wherein *Q* = *Qn Qs Qz* T and *q* = *qn qs qz* T are the concentrated load vector and the uniformly distributed load vector, respectively; *Qn*, *Qs* , and *Qz* are concentrated forces in the directions of the *n*-axis, *s*-axis, and *z*-axis, respectively; *qn*, *qs* , and *qz* 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).

$$\mathbf{W} = \begin{pmatrix} \mathbf{W}\_{nk} & \mathbf{W}\_{sk} & \mathbf{W}\_{\bar{z}k} \end{pmatrix}^{\mathrm{T}} = \mathbf{H}\_1 \mathbf{A} + \mathbf{H}\_2 \mathbf{A}' \tag{29}$$

where

$$\begin{aligned} \begin{bmatrix} \mathbf{H}\_{1} \end{bmatrix}\_{3 \times 10} &= \left\{ \begin{bmatrix} \mathbf{B}\_{1} \end{bmatrix}^{\mathrm{T}} \quad \begin{bmatrix} \mathbf{B}\_{2} \end{bmatrix}^{\mathrm{T}} \quad \begin{bmatrix} \mathbf{B}\_{3} \end{bmatrix}^{\mathrm{T}} \right\}^{\mathrm{T}} \end{aligned} \tag{30}$$

$$\begin{aligned} \begin{bmatrix} \mathbf{H}\_{2} \end{bmatrix}\_{3 \times 10} &= \left\{ \begin{array}{cc} \mathbf{0} \end{array} \begin{bmatrix} \mathbf{0} \end{bmatrix}\_{1 \times 10} ^ {\mathrm{T}} & \begin{bmatrix} \mathbf{0} \end{bmatrix}\_{1 \times 10} ^ {\mathrm{T}} & \begin{bmatrix} \mathbf{B}\_{4} \end{bmatrix}^ {\mathrm{T}} \right\}^{\mathrm{T}} \end{aligned} \tag{31}$$

$$\mathbf{B}\_1 = \begin{pmatrix} \cos \alpha & -\sin \alpha & 0 & -\rho\_n & 0 & -D\_n & \frac{a\_k \cos \alpha}{2} & 0 & 0 & 0 \\ & & & & & \ddots & \ddots \end{pmatrix} \tag{32}$$

$$\mathbf{B}\_2 = \begin{pmatrix} -\sin\alpha & -\cos\alpha & 0 & \rho\_\mathbf{s} & 0 & D\_\mathbf{s} & -\frac{a\_\mathbf{k}\sin\alpha}{2} & 0 & 0 & 0 & 0 \end{pmatrix} \tag{33}$$

$$\mathbf{B}\_2 = \begin{pmatrix} 0 & 0 & 1 - \mathbf{v}k\_0 & 0 & -\omega\_\mathbf{v} & 0 & 0 & \frac{a\_\mathbf{k}(1-\mathbf{x}k\_0)}{2} & \omega\_\mathbf{v} & \omega\_\mathbf{v} \end{pmatrix} \tag{34}$$

$$\mathbf{B}\_3 = \begin{pmatrix} 0 & 0 & 1 - \mathbf{x}k\_0 & 0 & -\omega & 0 & 0 & \frac{a\_k(1-\mathbf{x}k\_0)}{2} & \psi\_\mathbf{c} & \psi\_\mathbf{s} \end{pmatrix} \tag{34}$$

$$\mathbf{B}\_4 = \begin{pmatrix} -\omega & -(\omega + \omega k\_0) & 0 & 0 & 0 & -\omega \end{pmatrix} \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \end{pmatrix} \tag{35}$$

$$\mathbf{B}\_{4} = \begin{pmatrix} -\mathbf{x} & -(y + \omega k\_{0}) & 0 & 0 & 0 & -\omega\_{\rm d} & 0 & 0 & 0 & 0 \\ & & & & & & & \dots \end{pmatrix} \tag{35}$$

$$\mathbf{A} = \begin{pmatrix} \ \ u & v \ w & w & \theta & \beta & \theta\_{\mathrm{d}} & \mathfrak{Q}\_{\mathrm{x}} & \mathfrak{Q}\_{\mathrm{z}} f\_{\mathrm{c}} f\_{\mathrm{s}} \end{pmatrix}^{\mathrm{T}} \tag{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:

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

$$\mathbf{d}\_{i}^{\mathbf{e}} = \begin{pmatrix} u\_{i} & u\_{i} \ & v\_{i} & v\_{i} \ & w\_{i} & \theta\_{i} & \theta\_{i} & \theta\_{\mathrm{di}} & \theta\_{\mathrm{di}} & \mathbf{2}\Omega\_{\mathrm{zi}} & \mathbf{2}\Omega\_{\mathrm{zi}} & f\_{\mathrm{ci}} & f\_{\mathrm{si}} & \mathbf{) for } i = 1, 2 \end{pmatrix} \\ \text{for } i = 1, 2 \quad \text{(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:

$$\mathbf{K}d^{\mathbf{e}} = \mathbf{F} \tag{39}$$

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

$$\mathbf{K} = \int\_{\boldsymbol{l}\_{\rm c}} \mathbf{N}^{\rm T} \left[ \iint\_{A\_{\rm s}} \mathbf{B}\_{\rm s}^{\rm T} \mathbf{S}^{\rm T} \mathbf{E}\_{\rm s} \mathbf{S} \mathbf{B}\_{\rm s} \mathrm{d}a + \iint\_{A\_{\rm c}} \mathbf{B}\_{\rm c}^{\rm T} \mathbf{S}^{\rm T} \mathbf{E}\_{\rm c} \mathbf{S} \mathbf{B}\_{\rm c} \mathrm{d}a + \int\_{2\mathbf{b}\_{\rm s}} \mathbf{B}\_{\rm 1}^{\rm T} \boldsymbol{\rho}\_{\rm 1} \boldsymbol{\rho}\_{\rm 1} \mathbf{S}\_{\rm 1} \mathbf{B}\_{\rm 1} \mathrm{d}x + T\_{\rm R} \right] \mathrm{N} \mathrm{d}z \tag{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:

$$F = F\_{\mathbf{q}} + F\_{\xi} + F\_{\mathbf{sh}} \tag{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).

$$F\_{\mathbf{q}} = \int\_{l\_{\mathrm{e}}} \left( N\_{\mathrm{F}}^{\mathrm{T}} H\_{\mathrm{1}}^{\mathrm{T}} + N\_{\mathrm{F}}^{\mathrm{T}} H\_{\mathrm{2}}^{\mathrm{T}} \right) q \mathrm{d}z + \left( N\_{\mathrm{F}}^{\mathrm{T}} H\_{\mathrm{1}}^{\mathrm{T}} + N\_{\mathrm{F}}^{\mathrm{T}} H\_{\mathrm{2}}^{\mathrm{T}} \right) \mathrm{Q} \tag{42}$$

$$F\_{\xi} = -\int\_{l\_{\mathrm{c}}} \mathbf{N}^{T} \mathbf{B}\_{\mathrm{c}}^{T} \left( \iint\_{A\_{\mathrm{c}}} \mathbf{S}^{T} \boldsymbol{\xi}(t, t\_{0}) \boldsymbol{\sigma}\_{\mathrm{c}}(t\_{0}) \mathrm{d}a \right) \mathrm{d}z \tag{43}$$

$$F\_{\mathbf{sh}} = \int\_{l\_{\rm e}} \mathbf{N}^{T} \mathbf{B}\_{\mathbf{c}}^{T} \left( \iint\_{A\_{\rm c}} \mathbf{S}^{T} \mathbf{E}\_{\mathbf{ce}}(t, t\_{0}) \, \boldsymbol{\varepsilon}\_{\mathbf{c}}^{\mathbf{sh}}(t) \, \mathrm{d}a \right) \mathrm{d}z \tag{44}$$

where

$$\sigma\_{\mathbf{c}}(t\_0) = \begin{cases} \sigma\_{\mathbf{c}}(t\_0) & \tau\_{\mathbf{c}}(t\_0) \end{cases} \text{ \(} \mathbf{t}\_0\text{)}\text{ \)}^{\text{T}} \tag{45}$$

$$\varepsilon\_{\mathbf{c}}^{\mathbf{sh}}(t) = \begin{Bmatrix} \varepsilon\_{\mathbf{c}}^{\mathbf{sh}}(t) & 0 \end{Bmatrix}^{\mathrm{T}} \tag{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/mm3 ; and a vertical load of 10 kN/mm acts at the section centroid line.

FE model.

acts at the section centroid line.

where

w, are both equal to 100 N/mm3

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

d

**sh** *F N B S E* **sh c ce c** *ε*

Notably, the calculation methods for the characteristic geometrical parameters of the section

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

0

sh c *t t*

0 c 0 c 0 *t t t*

( ) ( ) *σ***c**

under torsion and distortion in the finite beam element model can be found in the literature [6,23].

 *z* **q F 1 F 2 F F 2** *F N H N H q N H N H Q*

T T T T T T T T

1

0 0 , d d

, d d

T

*t t t a z*

 *t t t a z*

0

T

(42)

(43)

(44)

(45)

(46)

; and a vertical load of 10 kN/mm

the concrete shrinkage effect, which can be calculated by Equation (44).

e c

*l A*

e c

**sh c** *ε*

*l A*

T T T

*F N B S* **c c**

T T T

e

*l*

**3. Numerical Validation of the Beam Element Model**

u and

**Figure 2.** Geometric dimensions and load layout of the curved composite box beam used for validation (unit: mm): (**a**) cross section, (**b**) elevation. **Figure 2.** Geometric dimensions and load layout of the curved composite box beam used for validation (unit: mm): (**a**) cross section, (**b**) elevation.

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 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 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. *Materials* **2020**, *12*, x FOR PEER REVIEW 10 of 17with Shell181 elements, whereas the shear connectors are modeled with Combin14 elements. The finite element model is shown in Figure 3.

**Figure 3.** Elaborate finite element model of a curved composite box beam established in ANSYS. **Figure 3.** Elaborate finite element model of a curved composite box beam established in ANSYS.

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 × 105 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*0 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 Figures 4–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 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 × 105 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*0 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 Figures 4–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.

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.

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. *Materials* **2020**, *12*, x FOR PEER REVIEW 11 of 17

*Materials* **2020**, *12*, x FOR PEER REVIEW 11 of 17

**Figure 4.** Deflection variation in curved composite box girder, with respect to time: (**a**) variation of mid-span deflection with respect to time, (**b**) distribution of deflection along the span. 0.3 7 **Figure 4.** Deflection variation in curved composite box girder, with respect to time: (**a**) variation of mid-span deflection with respect to time, (**b**) distribution of deflection along the span. **Figure 4.** Deflection variation in curved composite box girder, with respect to time: (**a**) variation of mid-span deflection with respect to time, (**b**) distribution of deflection along the span. (a) Variation of mid-span deflection with respect to time (b) Distribution of deflection along the span **Figure 4.** Deflection variation in curved composite box girder, with respect to time: (**a**) variation of mid-span deflection with respect to time, (**b**) distribution of deflection along the span.

(b) Distribution of deflection along the span

(a) Variation of mid-span deflection with respect to time

**Figure 5.** Transverse displacement variation in the curved composite box girder, with respect to (a) Variation of transverse mid-span displacement with respect to time (b) Distribution of transverse displacement along the span (a) Variation of transverse mid-span displacement with respect to time (b) Distribution of transverse displacement along the span

time: (**a**) variation of transverse mid-span displacement with respect to time, (**b**) distribution of transverse displacement along the span. 4 2 **Figure 5.** Transverse displacement variation in the curved composite box girder, with respect to time: (**a**) variation of transverse mid-span displacement with respect to time, (**b**) distribution of transverse displacement along the span. 4 2 **Figure 5.** Transverse displacement variation in the curved composite box girder, with respect to time: (**a**) variation of transverse mid-span displacement with respect to time, (**b**) distribution of transverse displacement along the span. **Figure 5.** Transverse displacement variation in the curved composite box girder, with respect to time: (**a**) variation of transverse mid-span displacement with respect to time, (**b**) distribution of transverse displacement along the span. 4 2

**Figure 6.** Stress variation in the steel girder and concrete slab of the curved composite box beam, with respect to time: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab. **Figure 6.** Stress variation in the steel girder and concrete slab of the curved composite box beam, with respect to time: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab. **Figure 6.** Stress variation in the steel girder and concrete slab of the curved composite box beam, with respect to time: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab. **Figure 6.** Stress variation in the steel girder and concrete slab of the curved composite box beam, with respect to time: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab.

**4. Parametric Analysis**

**4. Parametric Analysis**

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

**4. Parametric Analysis**

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.

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.

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.

#### **4. Parametric Analysis**

### *4.1. Influence of the Time-Dependent E*ff*ects 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 Figures 7 and 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. *Materials* **2020**, *12*, x FOR PEER REVIEW 12 of 17 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 Figures 7 and 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. *Materials* **2020**, *12*, x FOR PEER REVIEW 12 of 17 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 Figures 7 and 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.

**Figure 7.** Change in the shear-lag warping strength function of the concrete slab, with respect to time: (**a**) shear-lag coefficient at the inside; (**b**) shear-lag coefficient at the outside. **Figure 7.** Change in the shear-lag warping strength function of the concrete slab, with respect to time: (**a**) shear-lag coefficient at the inside; (**b**) shear-lag coefficient at the outside. **Figure 7.** Change in the shear-lag warping strength function of the concrete slab, with respect to time: (**a**) shear-lag coefficient at the inside; (**b**) shear-lag coefficient at the outside.

**Figure 8.** Change in the shear-lag warping strength function of the steel bottom flange, with respect to time: (**a**) shear-lag coefficient at the inside, (**b**) shear-lag coefficient at the outside. **Figure 8.** Change in the shear-lag warping strength function of the steel bottom flange, with respect to time: (**a**) shear-lag coefficient at the inside, (**b**) shear-lag coefficient at the outside. **Figure 8.** Change in the shear-lag warping strength function of the steel bottom flange, with respect to time: (**a**) shear-lag coefficient at the inside, (**b**) shear-lag coefficient at the outside.

#### *4.2. Influence of Time-Dependent Effects on Mechanical Behavior for Different Initial Curvatures 4.2. Influence of Time-Dependent Effects on Mechanical Behavior for Different Initial Curvatures 4.2. Influence of Time-Dependent E*ff*ects on Mechanical Behavior for Di*ff*erent 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 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 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

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.

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.

direction are calculated for the concrete slab and steel bottom flange at the initial and ultimate

and 60° are 58.8%, 57.8%, and 53.0%, respectively.

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. *Materials* **2020**, *12*, x FOR PEER REVIEW 13 of 17 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. *Materials* **2020**, *12*, x FOR PEER REVIEW 13 of 17 The compressive stress reduction rates at the inside of the concrete slab for central angles of 0°, 30°,

**Figure 9.** Displacement distribution along the span of the curved composite box girder, with different central angles at the initial loading age and ultimate loading age: (**a**) deflection along the span, (**b**)transverse displacement along the span, (**c**) longitudinal displacement along the span. 4 2 **Figure 9.** Displacement distribution along the span of the curved composite box girder, with different central angles at the initial loading age and ultimate loading age: (**a**) deflection along the span, (**b**)transverse displacement along the span, (**c**) longitudinal displacement along the span. **Figure 9.** Displacement distribution along the span of the curved composite box girder, with different central angles at the initial loading age and ultimate loading age: (**a**) deflection along the span, (**b**)transverse displacement along the span, (**c**) longitudinal displacement along the span.

girder, with different central angles at the initial loading age and ultimate loading age: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab. **Figure 10.** Normal stress distribution along the transverse direction of the curved composite box girder, with different central angles at the initial loading age and ultimate loading age: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab. **Figure 10.** Normal stress distribution along the transverse direction of the curved composite box girder, with different central angles at the initial loading age and ultimate loading age: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab.

#### 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. 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 *4.3. Influence of the Time-Dependent E*ff*ects on Mechanical Behavior for Di*ff*erent Interfacial Shear Connection Sti*ff*ness Values*

The compressive stress reduction rate of the concrete slab decreases gradually with increasing initial curvature. 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. 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/mm3 . According to the proposed finite beam element model, the deflection and transverse and longitudinal displacement distributions along

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

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

The curved composite box girder in Section 3 was selected as the research object. The interfacial

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

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

, respectively, due to the influence of concrete shrinkage and creep. Moreover,

, respectively, due to the influence of concrete shrinkage and creep. Moreover,

. According to the proposed finite

. According to the proposed finite

, the maximum growth rates of

, the maximum growth rates of

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

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

shear connection stiffness values are set as 0.1, 1, and 10 N/mm3

shear connection stiffness values are set as 0.1, 1, and 10 N/mm3

for shear connection stiffness values of 0.1, 1, and 10 N/mm3

for shear connection stiffness values of 0.1, 1, and 10 N/mm3

*Connection Stiffness Values*

*Connection Stiffness Values*

0.1, 1, and 10 N/mm3

0.1, 1, and 10 N/mm3

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/mm3 , respectively, due to the influence of concrete shrinkage and creep. Moreover, for shear connection stiffness values of 0.1, 1, and 10 N/mm3 , 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/mm3 , 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. *Materials* **2020**, *12*, x FOR PEER REVIEW 14 of 17 curved composite box girders with shear connection stiffness values of 0.1, 1, and 10 N/mm3 , 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. 1 2 *Materials* **2020**, *12*, x FOR PEER REVIEW 14 of 17 curved composite box girders with shear connection stiffness values of 0.1, 1, and 10 N/mm3 , 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%,

**Figure 11.** Displacement distribution along the span of curved composite box beam at the initial and ultimate loading ages for different interfacial shear connection stiffness: (**a**) deflection along the span, (**b**) Transverse displacement along the span, (**c**) longitudinal displacement along the span. 6 0 **Figure 11.** Displacement distribution along the span of curved composite box beam at the initial and ultimate loading ages for different interfacial shear connection stiffness: (**a**) deflection along the span, (**b**) Transverse displacement along the span, (**c**) longitudinal displacement along the span. **Figure 11.** Displacement distribution along the span of curved composite box beam at the initial and ultimate loading ages for different interfacial shear connection stiffness: (**a**) deflection along the span, (**b**) Transverse displacement along the span, (**c**) longitudinal displacement along the span.

girder at the initial and ultimate loading ages for different interfacial shear connection stiffness: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab. **Figure 12.** Normal stress distribution along the transverse direction of the curved composite box girder at the initial and ultimate loading ages for different interfacial shear connection stiffness: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab. **Figure 12.** Normal stress distribution along the transverse direction of the curved composite box girder at the initial and ultimate loading ages for different interfacial shear connection stiffness: (**a**) normal stress in the steel girder, (**b**) normal stress in the concrete slab.

#### 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 It can be concluded that the growth rate of displacement, which is caused by concrete **5. Conclusions**

**5. Conclusions**

**5. Conclusions**

elaborate finite element model.

elaborate finite element model.

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 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 (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

(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

(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

(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

(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

concrete slab becomes large with increasing interfacial shear connection stiffness.

the finite element method and the age-adjusted effective modulus method.

the finite element method and the age-adjusted effective modulus method.

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.

**Author Contributions:** Data curation, X.-L.J.; Formal analysis, G.-M.W. and L.Z.; Funding acquisition, L.Z.; Investigation, L.Z.; Resources, W.-Y.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Transportation Science and Technology Program of Hebei Province, China (grant number QG2018-8), the Science and Technology Program of Qinghai Province, China (grant number 2017-SF-139), and the National Natural Science Foundation of China (grant number 51608031).

**Acknowledgments:** The authors express their thanks to the people helping with this work, and acknowledge the valuable suggestions from the peer reviewers.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

A quadratic polynomial is used to characterize the shear-lag warping shape function, ψc(*x*), of the concrete slab and the shear-lag warping shape function, ψs(*x*), of the steel beam, which are expressed in Equations (A1) and (A2), respectively. The physical meaning of the polynomial is the shear-lag warping displacement at any position on the section of the concrete slab or steel beam when a unit shear-lag warping intensity occurs.

$$\psi\_{\mathbf{c}}(\mathbf{x}) = \begin{cases} \left[1 - \left(\frac{\mathbf{x} + b\_{\mathbf{c}}}{b\_{\mathbf{c}} - b\_{\mathbf{t}}}\right)^2\right] \left(\frac{b\_{\mathbf{c}} - b\_{\mathbf{t}}}{b\_{\mathbf{t}}}\right)^2 & \mathbf{x} > b\_{\mathbf{t}} \\\ 1 - \left(\frac{\mathbf{x}}{b\_{\mathbf{t}}}\right)^2 & -b\_{\mathbf{t}} \le \mathbf{x} \le b\_{\mathbf{t}} \\\ \left[1 - \left(\frac{\mathbf{x} - b\_{\mathbf{c}}}{b\_{\mathbf{c}} - b\_{\mathbf{t}}}\right)^2\right] \left(\frac{b\_{\mathbf{c}} - b\_{\mathbf{t}}}{b\_{\mathbf{t}}}\right)^2 & \mathbf{x} < -b\_{\mathbf{t}} \end{cases} \tag{A1}$$

$$\psi\_{\mathbf{s}}(\mathbf{x}) = \begin{cases} 0 & \text{top flame and web of the steel girder} \\ 1 - \left(\frac{\mathbf{x}}{b\_{\mathbf{s}}}\right)^2 & \text{bottom flame of the steel girder} \end{cases} \tag{A2}$$

### **Appendix B**

(1) The nonzero elements of matrix *Bk* are shown hereafter:

*Bk* (1,1) = −*k*0; *Bk* (1,8) = 1; *Bk* (1,9) = −*y*s*k*0; *Bk* (1,13) = −*y*d*k*0Ψd; *Bk* (1,16) = −*akk*0/2; *Bk* (1,19) = *ak* /2; *Bk* (2,3) = −1; *Bk* (2,8) = −*k*0; *Bk* (2,19) = −*akk*0/2;

*Bk* (3,6) = −1; *Bk* (3,9) = *k*0; *Bk* (3,13) = *k*0Ψd; *Bk* (4,6) = −*k*0; *Bk* (4,12) = −1; *Bk* (5,15) = −1; *Bk* (6,21) = 1; *Bk* (7,23) = 1; *Bk* (8,5) = *k*0; *Bk* (8,10) = 1; *Bk* (9,5) = −*k*0; *Bk* (9,11) = −1; *Bk* (10,20) = 1; and *Bk* (11,22) = 1. In these elements, *k* = s,c represents the steel girder and concrete slab, respectively. (2) The nonzero elements of matrix *B***Ω** are given hereafter: *B*Ω(1,16) = 1 and *B*Ω(2,18) = 1.

### **Appendix C**

The nonzero elements of the shape function matrix *N* are presented hereafter: *N*(1,1) = *n*1; *N*(1,2) = *n*2; *N*(1,14) = *n*3; *N*(1,15) = *n*4; *N*(2,1) = *n*1 0 ; *N*(2,2) = *n*2 0 ; *N*(2,14) = *n*3 0 ; *N*(2,15) = *n*4 0 ; *N*(3,1) = *n*1"; *N*(3,2) = *n*2"; *N*(3,14) = *n*3"; *N*(3,15) = *n*4"; *N*(4,3) = *n*1; *N*(4,4) = *n*2; *N*(4,16) = *n*3; *N*(4,17) = *n*4; *N*(5,3) = *n*1 0 ; *N*(5,4) = *n*2 0 ; *N*(5,16) = *n*3 0 ; *N*(5,17) = *n*4 0 ; *N*(6,3) = *n*1"; *N*(6,4) = *n*2"; *N*(6,16) = *n*3"; *N*(6,17) = *n*4"; *N*(7,5) = *m*1; *N*(7,18) = *m*2; *N*(8,5) = *m*1 0 ; *N*(8,18) = *m*2 0 ; *N*(9,6) = *m*1; *N*(9,19) = *m*2; *N*(10,6) = *m*1 0 ; *N*(10,19) = *m*2 0 ; *N*(11,7) = *m*1; *N*(11,20) = *m*2; *N*(12,7) = *m*1 0 ; *N*(12,20) = *m*2 0 ; *N*(13,8) = *n*1; *N*(13,9) = *n*2; *N*(13,21) = *n*3; *N*(13,22) = *n*4; *N*(14,8) = *n*1 0 ; *N*(14,9) = *n*2 0 ; *N*(14,21) = *n*3 0 ; *N*(14,22) = *n*4 0 ; *N*(15,8) = *n*1"; *N*(15,9) = *n*2"; *N*(15,21) = *n*3"; *N*(15,22) = *n*4"; *N*(16,10) = *m*1; *N*(16,23) = *m*2; *N*(17,10) = *m*1 0 ; *N*(17,23) = *m*2 0 ; *N*(18,11) = *m*1; *N*(18,24) = *m*2; *N*(19,11) = *m*1 0 ; *N*(19,24) = *m*2 0 ; *N*(20,12) = *m*1; *N*(20,25) = *m*2; *N*(21,12) = *m*1 0 ; *N*(21,25) = *m*2 0 ; *N*(22,13) = *m*1; *N*(22,26) = *m*2; *N*(23,13) = *m*1 0 ; and *N*(23,26) = *m*2 0 . The nonzero elements of the shape function matrix *N***F** are given hereafter: *N*F(1,1) = *n*1; *N*F(1,2) = *n*2; *N*F(1,14) = *n*3; *N*F(1,15) = *n*4; *N*F(2,3) = *n*1; *N*F(2,4) = *n*2; *N*F(2,16) = *n*3; *N*F(2,17) = *n*4; *N*F(3,5) = *m*1; *N*F(3,18) = *m*2; *N*F(4,6) = *m*1; *N*F(4,19) = *m*2; *N*F(5,7) = *m*1; *N*F(5,20) = *m*2; *N*F(6,8) = *n*1; *N*F(6,9) = *n*2; *N*F(6,21) = *n*3; *N*F(6,22) = *n*4; *N*F(7,10) = *m*1; *N*F(7,23) = *m*2; *N*F(8,11) = *m*1; *N*F(8,24) = *m*2; *N*F(9,12) = *m*1; *N*F(9,25) = *m*2; *N*F(10,13) = *m*1; and *N*F(10,26) = *m*2.

In these elements, *n*1 = 1 + 2*z l*e *z*−*l*e *l*e 2 ; *n*2 = *z z*−*l*e *l*e 2 ; *n*3 = 1 − 2(*z*−*l*e) *l*e *z l*e 2 ; *n*4 = (*z* − *l*e) *z l*e 2 ; *m*1 = 1 − *z l*e ; and *m*2 = *z l*e .

### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Axial Compressive Behaviour of Square Through-Beam Joints between CFST Columns and RC Beams with Multi-Layers of Steel Meshes**

#### **Weining Duan 1 , Jian Cai 1,2, Xu-Lin Tang 3 , Qing-Jun Chen 1,2, Chun Yang 1,\* and An He 4**

Received: 11 May 2020; Accepted: 27 May 2020; Published: 29 May 2020

**Abstract:** The axial compressive behaviour of an innovative type of square concrete filled steel tube (CFST) column to reinforced concrete (RC) beam joint was experimentally investigated in this paper. The innovative joint was designed such that (i) the steel tubes of the CFST columns were completely interrupted in the joint region, (ii) the longitudinal reinforcements from the RC beams could easily pass through the joint area and (iii) a reinforcement cage, including a series of reinforcement meshes and radial stirrups, was arranged in the joint area to strengthen the mechanical performance of the joint. A two-stage experimental study was conducted to investigate the behaviour of the innovative joint under axial compression loads, where the first stage of the tests included three full-scale innovative joint specimens subjected to axial compression to assess the feasibility of the joint detailing and propose measures to further improve its axial compressive behaviour, and the second stage of the tests involved 14 innovative joint specimens with the improved detailing to study the effect of the geometric size of the joint, concrete strength and volume ratio of the steel meshes on the bearing strengths of the joints. It was generally found from the experiments that (i) the innovative joint is capable of achieving the design criterion of the 'strong joint-weak member' with appropriate designs, and (ii) by decreasing the height factor and increasing the volume ratio of the steel meshes, the axial compressive strengths of the joints significantly increased, while the increase of the length factor is advantageous but limited to the resistances of the joint specimens. Because of the lack of existing design methods for the innovative joints, new design expressions were proposed to calculate the axial compression resistances of the innovative joints subjected to bearing loads, with the local compression effect, the confinement effect provided by the multi-layers of steel meshes and the height effect of concrete considered. It was found that the proposed design methods were capable of providing accurate and safe resistance predictions for the innovative joints.

**Keywords:** through-beam joint; concrete filled steel tube (CFST) columns; reinforced concrete (RC); axial compressive behaviour; steel mesh; local compression; confined concrete; height factor

## **1. Introduction**

The concrete-filled steel tube (CFST) structure is a high-efficiency solution for high-rise buildings and bridges due to its high bearing capacity, excellent ductility and great energy dissipation capacity [1].

In CFST structures, the joints connecting the CFST columns and reinforced concrete (RC) beams play an important role in achieving the structural integrity of the CFST frame [2,3], and thus researches related to the detailing of the joints have become a hot spot during the past few years. One of the typical joints connecting the CFST columns and RC beams is the through-column joint, in which steel cleats are attached to the outer surface of the steel tube of the CFST column at the joint area and the longitudinal reinforcements in the RC beams were welded on these steel cleats. The seismic behaviour and the failure modes of the through-column joints have been widely investigated [4–7], indicating that the through-column joints possess excellent seismic behaviour but require a significant amount of site welding.

Recently, through-beam joints connecting the CFST columns and RC beams have attracted the interest of researchers because of their simple detailing in the joint area, with a brief review of the recently developed through-beam joints summarised herein. Nie et al. [8] and Bai et al. [9] developed a representative through-beam joint connecting the concrete encased CFST columns and RC beams, where the steel tubes of the CFST columns were completely interrupted in the joint area, and the steel-reinforcement bars in the RC beams were continuous in the floor. Zhang et al. [10] developed a through-beam joint for concrete-filled double-skin steel tubes structural members, where the outer tube was interrupted in the joint area and an octagonal RC ring beam was arranged outside the joint to connect the reinforcements from the RC beams. Chen et al. [11,12] reported a new type of through-beam connection for CFST columns and RC beams, where the steel tube of the CFST column is completely or partially interrupted, and a ring beam was used to strengthen the load-carrying capacity of the joint. Zhou et al. [13,14] proposed a tubed-reinforced-concrete column to RC beam joint, where the steel tubes are completely interrupted in the joint area, with different manners to strengthen the joint, including the strengthening stirrup, horizontal hunches, and internal diaphragm. The seismic behaviour of the aforementioned through-beam joints was examined through cyclic loading tests. It was generally found that the through-beam joints are capable of achieving the design criterion of 'strong joint-weak member' with excellent seismic performance. Provided that the steel tubes of the CFST columns were generally interrupted at the joint area, axial compressive tests on the through-beam joints were also conducted by Nie et al. [8], Bai et al. [9], Chen et al. [11] and Zhou et al. [13], to investigate the axial compressive strengths of the joints; it was found that the strengths of the joints were higher than those of the corresponding CFST columns, and thus the design criterion of the 'strong joint-weak member' was deemed to be achieved. Recently, an innovative joint connecting the square CFST columns and RC beams has been developed by the authors [15]. The configuration and detailing of the innovative joint are displayed in Figure 1, where the outer steel tube was completely interrupted in the joint zone, and the longitudinal reinforcements from the RC beams could easily pass through the joint area. A reinforcement cage, including a series of reinforcement meshes and radial stirrups, was arranged in the joint area to strengthen the mechanical performance of the joint because of the interruption of the steel tube. The benefits of the proposed configuration of the innovative joints over the existing through-column and through-column beam joints include (i) the reinforcements from both structural columns and beams could easily pass through the joint area, ensuring both robustness and integrity of the joints, and (ii) the reinforcement cage in the joint can be precast and on-site welding is not required, leading to a simplified construction process. Cyclic loading tests on the innovative joints have been conducted by the authors [15] to study their seismic behaviour, while the axial compression behaviour of the innovative joints subjected to bearing loads remained unexplored and will be studied in the present paper.

Because of the lack of existing experimental investigation on the compressive behaviour of innovative joints, a two-stage experimental study was first conducted to investigate the compressive behaviour of the innovative joint subjected to bearing loads, where the first stage of the tests included three full-scale innovative joint specimens to assess the feasibility of the joint detailing and propose measures to further improve the axial compressive behaviour of the joint, and the second stage of the tests involved 14 innovative joint specimens with the improved detailing to further study the effects of

the geometric size of the joint, concrete strength and volume ratio of the steel meshes on the bearing strengths of the joints. Because of the lack of existing design methods for innovative joint, new design expressions for the bearing strengths of the innovative joints were developed and validated by the experimental results. *Materials* **2020**, *13*, x FOR PEER REVIEW 3 of 19 expressions for the bearing strengths of the innovative joints were developed and validated by the experimental results.

**Figure 1.** Square through-beam joint system.

#### **Figure 1.** Square through-beam joint system. **2. Tests of Joint with CFST Columns (Series I Tests)**

#### **2. Tests of Joint with CFST Columns (Series I Tests)**  *2.1. General*

SC2 500 × 6 1000 × 1000 ×

SC3 500 × 6 1000 × 1000 ×

(6)–(8).

*2.1. General*  The Series I axial compression tests included three geometrically identical full-scale specimens, namely SC1, SC2, and SC3, and were designed with various rebar diameters of the steel meshes in the joints. Each specimen included the upper and lower CFST columns and the joint, as shown in Figure 2, while the slab and the RC beams were not fabricated because of the size limitation of the testing machine. Note that the longitudinal reinforcements of the RC frame beams were still arranged in the joint. Figure 3 and Table 1 report the geometric dimensions and reinforcements of each part of the specimens. Specifically, the longitudinal reinforcements of the column were inserted through the joint area with the anchorage lengths of 1000 mm in the CFST columns, as determined based on the Chinese building code [16]. To avoid the compressive stress concentrating at the end of the steel tube, a steel ring whose width and thickness were 18 mm and 6 mm, respectively, was welded to each end of the steel tubes. Five layers of steel meshes with the adjacent distance of 75 mm were arranged in The Series I axial compression tests included three geometrically identical full-scale specimens, namely SC1, SC2, and SC3, and were designed with various rebar diameters of the steel meshes in the joints. Each specimen included the upper and lower CFST columns and the joint, as shown in Figure 2, while the slab and the RC beams were not fabricated because of the size limitation of the testing machine. Note that the longitudinal reinforcements of the RC frame beams were still arranged in the joint. Figure 3 and Table 1 report the geometric dimensions and reinforcements of each part of the specimens. Specifically, the longitudinal reinforcements of the column were inserted through the joint area with the anchorage lengths of 1000 mm in the CFST columns, as determined based on the Chinese building code [16]. To avoid the compressive stress concentrating at the end of the steel tube, a steel ring whose width and thickness were 18 mm and 6 mm, respectively, was welded to each end of the steel tubes. Five layers of steel meshes with the adjacent distance of 75 mm were arranged in the joint area. The diameters of the rebars in the steel meshes were 11 mm, 8.5 mm and 5.5 mm for specimens SC1, SC2 and SC3, respectively, as shown in Figure 3.

the joint area. The diameters of the rebars in the steel meshes were 11 mm, 8.5 mm and 5.5 mm for **Table 1.** Specimen parameters and the experimental and calculated loads for Series I tests.

**Specimen** *D* **×** *t* **(mm × mm)**  *A* **×** *A* **×** *H* **(mm × mm × mm)**  *f***cu,u (MPa)**  *f***cu (MPa)**  *f***cu,l (MPa)**  *d***bar (mm)**  ρ**v (%)**  *N***cr** **(kN)**  *N***u** **(kN)** *N***cr***/N***u** SC1 500 × 6 1000 × 1000 × 60046.3 40.8 44.7 11 1.65 6000 10,534 0.57 20,112 1.91 Notes: *D* and *t* are diameter and thickness of steel tube; *A* and *H* are width and height of the joint, respectively; *f* cu,u, *f* cu and *f* cu,l are the cubic compressive strengths of the concrete for upper column, joint area and lower column, respectively; *d*bar is the diameter of rebar in the steel meshes; ρv is the volume ratio of steel mesh in the joint, which is calculated as *V*sm/(*A* × *A* × *H*); *V*sm is the total volume of steel meshes in the joint; *N*cr is the load when the first crack appeared; *N*u is the peak load of the specimen; *N*cal is the calculated strength obtained by Equations (6)–(8).

60046.3 40.8 44.7 8.5 0.99 5300 12,004 0.44 18,395 1.53

60046.3 40.8 44.7 5.5 0.41 5100 11,404 0.45 15,954 1.40

*N***cal**

**(kN)** *N***cal***/N***u**

respectively; *d*bar is the diameter of rebar in the steel meshes; *ρ*v is the volume ratio of steel mesh in the joint, which is calculated as *V*sm/ (*A* × *A* × *H*); *V*sm is the total volume of steel meshes in the joint; *N*cr is the load when the first crack appeared; *N*u is the peak load of the specimen; *N*cal is the calculated strength obtained by Equations

*Materials* **2020**, *13*, x FOR PEER REVIEW 4 of 19

*Materials* **2020**, *13*, x FOR PEER REVIEW 4 of 19

**Figure 2.** Overview of the specimens. **Figure 2.** Overview of the specimens. **Figure 2.** Overview of the specimens.

**Figure 3.** Details of the specimens in Series I (**a**) Elevation view; (**b**) Top view; (**c**) Steel meshes and beam rebars in the joints (SG: strain gauge). **Figure 3.** Details of the specimens in Series I (**a**) Elevation view; (**b**) Top view; (**c**) Steel meshes and beam rebars in the joints (SG: strain gauge). **Figure 3.** Details of the specimens in Series I (**a**) Elevation view; (**b**) Top view; (**c**) Steel meshes andbeam rebars in the joints (SG: strain gauge).

During the fabrication of the specimens, four concrete cube specimens were respectively reserved when casting the concrete of the lower CFST column, joint area and upper CFST column. During the fabrication of the specimens, four concrete cube specimens were respectively reserved when casting the concrete of the lower CFST column, joint area and upper CFST column. During the fabrication of the specimens, four concrete cube specimens were respectively reserved when casting the concrete of the lower CFST column, joint area and upper CFST column. Compression

tests were conducted on the concrete cubes to derive their cubic strengths at the day of the column tests, with the average cubic compressive strength of each part of the specimen reported in Table 1. Tensile tests were conducted to derive the material properties of the reinforcements and the steel tube, with the key results from the tensile tests reported in Table 2. The Young's modulus of the steels may be taken as 200,000 MPa, according to the Chinese building code [16]. *Materials* **2020**, *13*, x FOR PEER REVIEW 5 of 19 in Table 1. Tensile tests were conducted to derive the material properties of the reinforcements and

**Table 2.** Measured material properties of reinforcing bars and steel tube for Series I test. the steel tube, with the key results from the tensile tests reported in Table 2. The Young's modulus of the steels may be taken as 200,000 MPa, according to the Chinese building code [16].

#### *2.2. Experimental Test Setup* Ultimate strength *f*u (MPa) 578 460 608 540 575 575 460

its maximum strength.

The Changchun CSS-254 15 MN universal testing machine was employed to apply axial compression force to the specimens, with the test setup shown in Figure 4. During the tests, four LVDTs were arranged in the joint area and the CFST columns to measured axial shortening of the specimen (see Figure 4), and a series of strain gauges were attached to the steel meshes in the joint to measure their deformations during loading (See Figure 3). Both load and displacement control modes were employed to drive the testing machine. Specifically, the initial loading rate was set to be equal to 5 kN/s until the material yielding was observed, according to the readings from the attached strain gauges, after which a displacement rate of 2 mm/min was used to complete the post-yield loading stage. The tests were terminated when the resistance of the specimen had deteriorated below 85% of its maximum strength. *2.2. Experimental Test Setup*  The Changchun CSS-254 15 MN universal testing machine was employed to apply axial compression force to the specimens, with the test setup shown in Figure 4. During the tests, four LVDTs were arranged in the joint area and the CFST columns to measured axial shortening of the specimen (see Figure 4), and a series of strain gauges were attached to the steel meshes in the joint to measure their deformations during loading (See Figure 3). Both load and displacement control modes were employed to drive the testing machine. Specifically, the initial loading rate was set to be equal to 5 kN/s until the material yielding was observed, according to the readings from the attached strain gauges, after which a displacement rate of 2 mm/min was used to complete the post-yield loading stage. The tests were terminated when the resistance of the specimen had deteriorated below 85% of

**Figure 4.** Schematic diagram of the test setup. **Figure 4.** Schematic diagram of the test setup.

#### *2.3. Experimental Results 2.3. Experimental Results*

#### 2.3.1. Failure Modes 2.3.1. Failure Modes

All the specimens exhibited similar overall behaviour under axial compression during the tests, as summarised below. The first vertical crack was found on the side surface of the joint when the upper load reached 44~57% of the peak loads. With the load increased, the horizontal cracks occurred and gradually crossed the first vertical crack. The dilation and local buckling of the steel tubes were observed when the specimens reached their peak loads, while the concrete in the joint area remained All the specimens exhibited similar overall behaviour under axial compression during the tests, as summarised below. The first vertical crack was found on the side surface of the joint when the upper load reached 44~57% of the peak loads. With the load increased, the horizontal cracks occurred and gradually crossed the first vertical crack. The dilation and local buckling of the steel tubes were observed when the specimens reached their peak loads, while the concrete in the joint area remained

member'. The crack patterns on each side surface of the joint at failure are shown in Figure 5.

intact without series spalling, which indicates that the specimens were failed by the local buckling of outer steel tubes of the CFST columns and thus satisfied the design criterion of the 'strong joint-weak member'. The crack patterns on each side surface of the joint at failure are shown in Figure *Materials*  5. **2020**, *13*, x FOR PEER REVIEW 6 of 19

Surface-2 Surface-1 Surface-4 Surface-3

**Figure 5.** Cracks on the side surfaces of the specimens. **Figure 5.** Cracks on the side surfaces of the specimens.

Taken specimen SC2 as an example to illustrate the experimental observation in detail, when the load reached about 45% of *N*u (*N*u is the peak load of the specimen, for SC2 *N*u = 12,004 kN), the first crack was found at the middle of the side surface of the joint in the vertical direction. As the load increased to about 50% of *N*u, the vertical crack developed with a maximum width of about 0.1 mm, and a short horizontal crack occurred on the side surface of the concrete at the same time. At about 70% of *N*u, the maximum widths of the vertical and horizontal cracks on the side surface of the joint reached 0.25 mm and 0.15 mm, respectively. When the specimen reached its maximum load *N*u, dilation and local buckling of the steel tubes of the CFST columns were observed. At the post-peak loading stage, the steel tubes of the CFST columns continued to dilate with the load capacity of the specimen gradually dropped to about 85% *N*u. Taken specimen SC2 as an example to illustrate the experimental observation in detail, when the load reached about 45% of *N*u (*N*u is the peak load of the specimen, for SC2 *N*u = 12,004 kN), the first crack was found at the middle of the side surface of the joint in the vertical direction. As the load increased to about 50% of *N*u, the vertical crack developed with a maximum width of about 0.1 mm, and a short horizontal crack occurred on the side surface of the concrete at the same time. At about 70% of *N*u, the maximum widths of the vertical and horizontal cracks on the side surface of the joint reached 0.25 mm and 0.15 mm, respectively. When the specimen reached its maximum load *N*u, dilation and local buckling of the steel tubes of the CFST columns were observed. At the post-peak loading stage, the steel tubes of the CFST columns continued to dilate with the load capacity of the specimen gradually dropped to about 85% *N*u.

#### 2.3.2. Load-Deformation and Load-Strain Curves 2.3.2. Load-Deformation and Load-Strain Curves

The axial load-deformation curves of the three specimens are shown in Figure 6, where the vertical axis is the force applied at two ends of the columns, and the horizontal axis is the displacement of the whole specimen consisted of the upper and lower CFST columns and the joint zone. As shown in Figure 6, the initial elastic stiffnesses and the shapes of the curves for the three specimens are quite close to each other at the elastic stages. After reaching the peak loads of the specimens, the load-deformation curves descend gradually and show ductile post-peak behaviour. The loads corresponding to the first cracks *N*cr and the peak loads *N*u for all specimens are listed in Table 1. It can be seen that the value of *N*cr for SC1 is 1.13 and 1.18 times than those of SC2 and SC3 respectively, which indicates that the crack resistance of the joint is enhanced by improving the volume ratio of steel meshes in the joint (*ρ*v). However, the ultimate axial compressive strength of SC1 is smaller than those of SC2 and SC3. It may be attributed to the discreteness of concrete material in the CFST columns, as all specimens eventually failed in the CFST columns. The axial load-deformation curves of the three specimens are shown in Figure 6, where the vertical axis is the force applied at two ends of the columns, and the horizontal axis is the displacement of the whole specimen consisted of the upper and lower CFST columns and the joint zone. As shown in Figure 6, the initial elastic stiffnesses and the shapes of the curves for the three specimens are quite close to each other at the elastic stages. After reaching the peak loads of the specimens, the load-deformation curves descend gradually and show ductile post-peak behaviour. The loads corresponding to the first cracks *N*cr and the peak loads *N*u for all specimens are listed in Table 1. It can be seen that the value of *N*cr for SC1 is 1.13 and 1.18 times than those of SC2 and SC3 respectively, which indicates that the crack resistance of the joint is enhanced by improving the volume ratio of steel meshes in the joint (ρv). However, the ultimate axial compressive strength of SC1 is smaller than those of SC2 and SC3. It may be attributed to the discreteness of concrete material in the CFST columns, as all specimens eventually failed in the CFST columns.

The load-strain relationships of the steel meshes in all specimens are similar, as Figure 7 shows. It can be found that the strains of the steel meshes grew slowly in the early loading phase. As the upper load increased to about 50% of *N*u when the first crack appeared, the strains developed dramatically and subsequently reached the yield strains at the peak loads. It was also found that as the volume ratio of the steel meshes in the joint decreases, the strain values corresponding to the peak loads increase, which exhibits more adverse elongation deformation of the steel meshes.

*Materials* **2020**, *13*, x FOR PEER REVIEW 7 of 19