**4. Establishment of Numerical Modeling**

### *4.1. Modeling*

The specimen had two spans of 3 m and the total length of the girder was 6.3 m. A 3D numerical model of the steel–concrete composite girder beam with a clear span of 3 m was simulated in ABAQUS version 6.14 [50]. A detailed view of the numerical model is shown in Figure 6. Loading was applied to the top surface of the concrete slab and distributed over the full width of the girder. The load increased with 20 kN intervals up to 100 kN and, then, the load interval was set to 50 kN up to the total imposing force of the actuator. Due to the symmetry in geometry, loading, and boundary conditions, only half of the beam was modeled. The coordinate system represented axes X, Y, and Z as axes 1, 2, and 3 in the model, respectively. The symmetry boundary conditions are shown in Figure 6e with a restrained degree of freedom [51]. **4. Establishment of Numerical Modeling**  *4.1. Modeling*  The specimen had two spans of 3 m and the total length of the girder was 6.3 m. A 3D numerical model of the steel–concrete composite girder beam with a clear span of 3 m was simulated in ABAQUS version 6.14 [50]. A detailed view of the numerical model is shown in Figure 6. Loading was applied to the top surface of the concrete slab and distributed over the full width of the girder. The load increased with 20 kN intervals up to 100 kN and, then, the load interval was set to 50 kN up to the total imposing force of the actuator. Due to the symmetry in geometry, loading, and boundary

*Materials* **2019**, *12*, x FOR PEER REVIEW 10 of 26

ABAQUS/standard solver with a linear geometric order was used in the study. Eight-node brick elements with reduced integration (C3D8R) were used to model the concrete and stud. The first-order interpolation, three-dimensional beam element (B31), and four-node shell element with reduced integration (S4R) were used to model reinforcement bars and steel, respectively. The final mesh included 55,496 elements and 67,327 nodes. Steel reinforcement bars were modeled using embedded rebar element. Surface-to-surface contact was used for stud-to-concrete and steel-to-concrete interactions. The material nonlinearities of concrete and steel were modeled using the concrete damage plasticity model and the elastic-plastic bilinear model, respectively. The accuracy of the results basically depends on the size of the mesh, the constitutive model, and boundary conditions [52]. Therefore, these aspects should be carefully incorporated into the proposed finite element model. Adequate attention was paid to the development of hexahedral mesh and the assigning interaction between various surfaces. Various components, namely, concrete slab, steel beam, stud connectors, reinforcement bars, and stiffeners, were meshed part by part instead of using global or sweep features. Thus, a regular structured hexahedral mesh was generated. To get an acceptable level of accuracy, the approximate global mesh size of 0.025 was used for reinforcement bars and concrete; whereas an approximate global mesh size of 0.015 was used for studs and steels. The modeling of different parts of the beam is shown in Figure 6. Material and geometrical details of the specimen girder are given in Table 1. conditions, only half of the beam was modeled. The coordinate system represented axes X, Y, and Z as axes 1, 2, and 3 in the model, respectively. The symmetry boundary conditions are shown in Figure 6e with a restrained degree of freedom [51]. ABAQUS/standard solver with a linear geometric order was used in the study. Eight-node brick elements with reduced integration (C3D8R) were used to model the concrete and stud. The first-order interpolation, three-dimensional beam element (B31), and four-node shell element with reduced integration (S4R) were used to model reinforcement bars and steel, respectively. The final mesh included 55,496 elements and 67,327 nodes. Steel reinforcement bars were modeled using embedded rebar element. Surface-to-surface contact was used for stud-to-concrete and steel-to-concrete interactions. The material nonlinearities of concrete and steel were modeled using the concrete damage plasticity model and the elastic-plastic bilinear model, respectively. The accuracy of the results basically depends on the size of the mesh, the constitutive model, and boundary conditions [52]. Therefore, these aspects should be carefully incorporated into the proposed finite element model. Adequate attention was paid to the development of hexahedral mesh and the assigning interaction between various surfaces. Various components, namely, concrete slab, steel beam, stud connectors, reinforcement bars, and stiffeners, were meshed part by part instead of using global or sweep features. Thus, a regular structured hexahedral mesh was generated. To get an acceptable level of accuracy, the approximate global mesh size of 0.025 was used for reinforcement bars and concrete; whereas an approximate global mesh size of 0.015 was used for studs and steels. The modeling of different parts of the beam is shown in Figure 6. Material and geometrical details of the specimen girder are given in Table 1.

*Materials* **2019**, *12*, x FOR PEER REVIEW 11 of 26

**Figure 6.** The finite element model: (**a**) Reinforced-steel meshing; (**b**) concrete-slab meshing; (**c**) steel beam, diaphragm, and stud meshing; (**d**) composite box girder meshing; (**e**) boundary conditions; (**f**) load-deflection behavior at 1000 kN load. **Figure 6.** The finite element model: (**a**) Reinforced-steel meshing; (**b**) concrete-slab meshing; (**c**) steel beam, diaphragm, and stud meshing; (**d**) composite box girder meshing; (**e**) boundary conditions; (**f**) load-deflection behavior at 1000 kN load.

#### *4.2. Material Description 4.2. Material Description*

#### 4.2.1. Concrete 4.2.1. Concrete

In order to simulate the mechanical behavior of concrete, a damage plastic model of concrete was utilized in the current study. The constitutive relationship of concrete is adopted from the uniaxial tension and compressive stress–strain curve of concrete given in the "code for the design of concrete structure" GB 50010-2010 [53]. On the day of the experiment, the measured average compressive strength of C60 concrete cube was 76.24 MPa, and the design strength of reserved-hole C80 concrete was 87.4 MPa. Additionally, the average splitting strength and the average elastic modulus were 5.09 MPa and 36.4 GPa, respectively. The uniaxial tension and the uniaxial compression parameters are listed in Table 4. The stress–strain relationship of concrete under uniaxial tension is expressed as follows: In order to simulate the mechanical behavior of concrete, a damage plastic model of concrete was utilized in the current study. The constitutive relationship of concrete is adopted from the uniaxial tension and compressive stress–strain curve of concrete given in the "code for the design of concrete structure" GB 50010-2010 [53]. On the day of the experiment, the measured average compressive strength of C60 concrete cube was 76.24 MPa, and the design strength of reserved-hole C80 concrete was 87.4 MPa. Additionally, the average splitting strength and the average elastic modulus were 5.09 MPa and 36.4 GPa, respectively. The uniaxial tension and the uniaxial compression parameters are listed in Table 4.

 = ( 1 − ௧) (1) **Table 4.** Concrete uniaxial tension and compression parameters.


The stress–strain relationship of concrete under uniaxial tension is expressed as follows:

$$
\sigma = (1 - d\_t) E\_\varepsilon \varepsilon \tag{1}
$$

$$d\_l = 1 - \rho\_l \left[ 1.2 - 0.2x^5 \right]\_{l'} \ge 1 \tag{2}$$

$$d\_l = 1 - \frac{\rho\_l}{a\_l(\mathbf{x} - \mathbf{1})^{1.7} + \mathbf{x}}, \mathbf{x} > 1\tag{3}$$

$$\mathfrak{x} = \frac{\mathfrak{e}}{\mathfrak{e}\_{t,r}} \tag{4}$$

$$
\rho\_l = \frac{f\_{l,r}}{E\_c \varepsilon\_{l,r}} \tag{5}
$$

where *α<sup>t</sup>* is the parameter of concrete uniaxial tension stress–strain curve in the decline period, *ft*,*<sup>r</sup>* is the representation of concrete uniaxial tensile strength, *εt*,*<sup>r</sup>* is the peak tensile strain corresponding to *ft*,*r* , and *d<sup>t</sup>* is the evolution parameter of concrete under uniaxial tension. to ௧,, and ௧ is the evolution parameter of concrete under uniaxial tension. Compression stress–strain relationships are assumed as follows: where is the parameter of concrete uniaxial compression stress–strain curve in the decline period,

Compression stress–strain relationships are assumed as follows: = ( 1 − ) (6)

$$
\sigma = (1 - d\_{\mathfrak{c}}) E\_{\mathfrak{c}} \mathfrak{e} \tag{6}
$$

$$d\_c = 1 - \frac{\rho\_c n}{n - 1 + \mathbf{x}^n}, \mathbf{x} \le 1 \tag{7}$$

$$\rho\_c$$

$$d\_c = 1 - \frac{\rho\_c}{a\_c(x-1)^2 + x}, x > 1\tag{8}$$

$$f\_{c,r}$$

$$\rho\_c = \frac{f\_{c,r}}{E\_c \varepsilon\_{c,r}} \tag{9}$$

$$n = \frac{E\_{\mathcal{C}} \varepsilon\_{\mathcal{C},r}}{E\_{\mathcal{C}} \varepsilon\_{\mathcal{C},r} - f\_{\mathcal{C},r}} \tag{10}$$

$$\mathbf{x} = \frac{\varepsilon}{\varepsilon\_{\mathcal{C},r}} \tag{11}$$

where *α<sup>c</sup>* is the parameter of concrete uniaxial compression stress–strain curve in the decline period, *fc*,*<sup>r</sup>* is the representation of concrete uniaxial compressive strength, *εt*,*<sup>r</sup>* is the peak compressive strain corresponding to *fc*,*<sup>r</sup>* , and *d<sup>c</sup>* is the evolution parameter of concrete under uniaxial compression. , is the representation of concrete uniaxial compressive strength, ௧, is the peak compressive strain corresponding to , , and is the evolution parameter of concrete under uniaxial compression. The stress–strain relationship and damage model of C60 concrete are shown in Figure 7.

The stress–strain relationship and damage model of C60 concrete are shown in Figure 7.

**Figure 7.** The stress–strain curve of C60 Concrete. (**a**) Compression; (**b**) tension. **Figure 7.** The stress–strain curve of C60 Concrete. (**a**) Compression; (**b**) tension.

**Table 4.** Concrete uniaxial tension and compression parameters.

elastic-plastic material model was employed based on the nominal stress–strain behavior of steel. In the stud connection, an elastic-plastic bilinear model was utilized. However, the material of steel

#### 4.2.2. Steel Beam, Stiffener, and Stud

**Parameter**  , (**Mpa**) , , (**Mpa**) , C60 3.82 3.5 1.28 × 10-4 3.81 76.24 2.19 × 10−<sup>3</sup> 4.2.2. Steel Beam, Stiffener, and Stud The steel and stiffener were modeled considering the nonlinear behavior of the materials. The The steel and stiffener were modeled considering the nonlinear behavior of the materials. The elastic-plastic material model was employed based on the nominal stress–strain behavior of steel. In the stud connection, an elastic-plastic bilinear model was utilized. However, the material of steel plates, studs, and steel bars was defined by the ideal elastoplastic model. That is, when the steel yields, the bearing capacity does not increase, but the deformation continues to increase. The stress–strain relationships of steel plates, studs, and bars are shown in Figure 8.

plates, studs, and steel bars was defined by the ideal elastoplastic model. That is, when the steel

**Figure 8.** Stress–strain relationships: (**a**) The stress–strain curve of steel; (**b**) the stress–strain curve of stud; (**c**) the stress–strain curve of ɸ8 rebar; (**d**) the stress–strain curve of ɸ10 rebar. **Figure 8.** Stress–strain relationships: (**a**) The stress–strain curve of steel; (**b**) the stress–strain curve of stud; (**c**) the stress–strain curve of φ8 rebar; (**d**) the stress–strain curve of φ10 rebar.

#### *4.3. Bonding 4.3. Bonding*

The bonding between the materials was done by the use of interaction in Abaqus. The stud and concrete was modeled by the penalty method considering a friction coefficient of 0.4 in the tangential direction and hard contact in the normal direction to avoid penetration between the two contact surfaces [51]. On account of the interaction between the flange of steel and concrete slabs, the steel was determined as the "slave surface" and the concrete as the "master surface". The finite sliding method was employed for the interaction between studs and concrete. To simulate the steel bar– concrete interaction, the reinforcement bar was selected as the embedded region and concrete was set to be the host region. The bonding between the materials was done by the use of interaction in Abaqus. The stud and concrete was modeled by the penalty method considering a friction coefficient of 0.4 in the tangential direction and hard contact in the normal direction to avoid penetration between the two contact surfaces [51]. On account of the interaction between the flange of steel and concrete slabs, the steel was determined as the "slave surface" and the concrete as the "master surface". The finite sliding method was employed for the interaction between studs and concrete. To simulate the steel bar–concrete interaction, the reinforcement bar was selected as the embedded region and concrete was set to be the host region.

#### *4.4. Comparison Between Numerical Analysis Values and Experimental Results 4.4. Comparison Between Numerical Analysis Values and Experimental Results*

The validation of experimental results was performed using the numerical analysis data as described earlier in Section 4. There was a good agreement between numerical and experimental results. Due to the maximum capacity of instrument available in the lab, a set of two actuators with a total imposing load of 900 kN was applied to the model. However, in the numerical model, the final step of the model was set to 1000 kN with the full shear interaction. The load increment was considered 10% in each load step. The numerical results showed a relatively stiffer behavior compared to the experimental model. It could be due to the stabilization process that occurred during the loading stage because of the change in loading system from one set of actuators to another set of actuators. However, similar results were obtained at the 900 kN load. Afterward, by a small The validation of experimental results was performed using the numerical analysis data as described earlier in Section 4. There was a good agreement between numerical and experimental results. Due to the maximum capacity of instrument available in the lab, a set of two actuators with a total imposing load of 900 kN was applied to the model. However, in the numerical model, the final step of the model was set to 1000 kN with the full shear interaction. The load increment was considered 10% in each load step. The numerical results showed a relatively stiffer behavior compared to the experimental model. It could be due to the stabilization process that occurred during the loading stage because of the change in loading system from one set of actuators to another set of actuators. However, similar results were obtained at the 900 kN load. Afterward, by a small increment in load, a significant change was observed in the deflection that resulted from the numerical model. Due to the symmetry in the geometric shape of the beam, the experimental measurement was performed in just one span.

The deflection at 900 kN was analytically found as 12.82 mm, which is very close to the experimental value of 13.211 mm. The deflection obtained from the numerical analysis is presented in Table 5. 12.82 mm, which is very close to the experimental value of 13.211 mm. The deflection obtained from the numerical analysis is presented in Table 5.

*Materials* **2019**, *12*, x FOR PEER REVIEW 14 of 26

increment in load, a significant change was observed in the deflection that resulted from the


**Table 5.** Load-deflection comparison of beam. **Table 5.** Load-deflection comparison of beam.

#### **5. Simplified Model 1** (**mm**)

The first crack observed when the load exceeded 80 kN; then, by increasing the load intensity, further crack formations were observed in the specimen. The cracks formed on the box girder were marked by a set of numbers demonstrating the occurrence order. Crack number 3 had the maximum crack width and was considered in the formulation process. The major cracks observed in the specimen are presented in Figures 9 and 10. Additionally, the related data are listed in Tables 2 and 3. According to EC4 [49], the effective length of the beam was considered 1500 mm. To observe the crack width and length and measure them accurately, a 5 × 5 cm grid size was made on concrete slab as demonstrated in Figures 9 and 10. The relationship between crack width and maximum central deflection of the ACHPCBG-bridge is addressed in this section. This relationship relies on the evaluation of load-deflection behaviors and load-crack width behaviors of experimental model outcomes. **5. Simplified Model**  The first crack observed when the load exceeded 80 kN; then, by increasing the load intensity, further crack formations were observed in the specimen. The cracks formed on the box girder were marked by a set of numbers demonstrating the occurrence order. Crack number 3 had the maximum crack width and was considered in the formulation process. The major cracks observed in the specimen are presented in Figures 9 and 10. Additionally, the related data are listed in Tables 2 and 3. According to EC4 [49], the effective length of the beam was considered 1500 mm. To observe the crack width and length and measure them accurately, a 5 × 5 cm grid size was made on concrete slab as demonstrated in Figures 9 and 10. The relationship between crack width and maximum central

There are many situations where an individual wants to use a simplified model and finds a formula that best fits a given set of data. Simplifying the model is the best solution and is the process of constructing a mathematical function with the best fit to a series of data points, giving a mathematical ideal solution. deflection of the ACHPCBG-bridge is addressed in this section. This relationship relies on the evaluation of load-deflection behaviors and load-crack width behaviors of experimental model outcomes.

(**b**)

**Figure 9.** Formation and distribution of cracks in specimen 1: (**a**) Schematic view of cracks; (**b**) actual view of cracks. **Figure 9.** Formation and distribution of cracks in specimen 1: (**a**) Schematic view of cracks; (**b**) actual view of cracks.

mathematical ideal solution.

(**b**)

**Figure 10.** Formation and distribution of cracks in specimen 2: (**a**) Schematic view of cracks; (**b**) actual view of cracks. **Figure 10.** Formation and distribution of cracks in specimen 2: (**a**) Schematic view of cracks; (**b**) actual view of cracks.

There are many situations where an individual wants to use a simplified model and finds a formula that best fits a given set of data. Simplifying the model is the best solution and is the process of constructing a mathematical function with the best fit to a series of data points, giving a As the first step, with experimental results a simplified formula was developed between the deflection and the load applied on the beam. It was done by fitting a first-order polynomial function to the data presented in Figure 11. The result was as follows:

$$L = 72.183d - 36.781\tag{12}$$

deflection and the load applied on the beam. It was done by fitting a first-order polynomial function to the data presented in Figure 11. The result was as follows: where (m) is the deflection measured at upper yield point of elastic stage in the center of the beam where *d* (m) is the deflection measured at upper yield point of elastic stage in the center of the beam and *L* (kN) is the intensity of the load excreted on the beam. *Materials* **2019**, *12*, x FOR PEER REVIEW 16 of 26

**Figure 11.** Fitting the load–deflection curve. **Figure 11.** Fitting the load–deflection curve.

Experimental results Proposed equation

**Figure 12.** Fitting the load–crack width curve.

<sup>0</sup> <sup>100</sup> <sup>200</sup> <sup>300</sup> <sup>400</sup> <sup>500</sup> <sup>600</sup> <sup>700</sup> <sup>800</sup> <sup>900</sup> <sup>1000</sup> <sup>0</sup>

Load (kN)

moment region is shown in Figure 13.

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

Crack width (mm)

To adopt the proposed model for the prediction of crack size in steel–concrete composite small box girder, one can simply calculate the maximum deflection under the static load case**.** The equivalent intensity of the load can then be calculated using Equation (12). Next, Equation (13) can be utilized to approximate the maximum width of the crack generated in box girder beam. However, the limitation of the proposed model regarding the dimension and construction size should also be considered while utilizing the proposed model. Having Equation (14), the crack width is available at any desired beam deflection. The behavior of central deflection and crack width at negative bending

**6. Discussion** 

100 200

Experimental results Proposed equation

Next, using the data obtained from the experimental tests, a third-order polynomial regression model was developed between the load and crack width. It was done by fitting a third-order polynomial function to the data presented in Figure 12. The result was as follows: 300 400 500 Load (kN)

*Materials* **2019**, *12*, x FOR PEER REVIEW 16 of 26

$$\mathcal{C} = 4 \ast 10^{-10} L^3 - 9 \ast 10^{-7} L^2 + 0.001 L - 0.0317 \tag{13}$$

where *C* (mm) is the maximum width of crack observed before failure. Substituting Equation (12) in Equation (13) gives a new expression, as follows: <sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> <sup>12</sup> <sup>14</sup> <sup>0</sup> Deflection (mm)

$$\mathcal{C} = 1.5 \ast 10^{-4} d^3 - 4.9193 \ast 10^{-3} d^2 + 0.07707 d - 0.06971 \tag{14}$$

**Figure 12.** Fitting the load–crack width curve. **Figure 12.** Fitting the load–crack width curve.

To adopt the proposed model for the prediction of crack size in steel–concrete composite small box girder, one can simply calculate the maximum deflection under the static load case**.** The equivalent intensity of the load can then be calculated using Equation (12). Next, Equation (13) can be utilized to approximate the maximum width of the crack generated in box girder beam. However, the limitation of the proposed model regarding the dimension and construction size should also be considered while utilizing the proposed model. Having Equation (14), the crack width is available at any desired beam deflection. The behavior of central deflection and crack width at negative bending moment region is shown in Figure 13. To adopt the proposed model for the prediction of crack size in steel–concrete composite small box girder, one can simply calculate the maximum deflection under the static load case. The equivalent intensity of the load can then be calculated using Equation (12). Next, Equation (13) can be utilized to approximate the maximum width of the crack generated in box girder beam. However, the limitation of the proposed model regarding the dimension and construction size should also be considered while utilizing the proposed model. Having Equation (14), the crack width is available at any desired beam deflection. The behavior of central deflection and crack width at negative bending moment region is shown in Figure 13. *Materials* **2019**, *12*, x FOR PEER REVIEW 17 of 26

**Figure 13.** Crack width-deflection comparison of beam. **Figure 13.** Crack width-deflection comparison of beam.

model was created. The experimental model was performed under different loads, and the maximum deflection was calculated in each case. Additionally, the crack width was measured from the experimental test. The results are reported in Tables 6 and 7. The load-deflation comparisons of tested specimen are shown in Figure 14. Furthermore, the crack width was approximated by the proposed model for each load intensity and compared with the experimental results of specimen 1 and 2.

**Table 6.** Validation of the proposed relationship for crack width prediction of specimen 1.

80 1.6 0.045 0.041 150 2.456 0.08 0.092 300 4.535 0.215 0.192 450 7.102 0.336 0.283 600 8.012 0.389 0.31 800 11.71 0.403 0.399 900 13.211 0.456 0.435

**Table 7.** Validation of the proposed relationship for crack width prediction of specimen 2.

80 1.2975 0.054 0.0223 150 2.456 0.094 0.0921 300 4.741 0.188 0.201 450 7.177 0.223 0.285 600 9.1455 0.265 0.338 800 11.8885 0.336 0.403 900 13.0965 0.349 0.432

**Crack width, mm**  (**Experimental model**)

**Crack Width, mm**  (**Experimental model**)

**Crack width, mm**  (**Proposed formula**)

> **Crack width, mm**  (**Proposed formula**)

**Load, KN Deflection, mm** 

**Load, KN Deflection, mm** 

**Experiment** 

**Experiment** 
