Analysis of Residual Flexural Stiffness of Steel Fiber-Reinforced Concrete Beams with Steel Reinforcement

This paper investigates the ability of steel fibers to enhance the short-term behavior and flexural performance of realistic steel fiber-reinforced concrete (SFRC) structural members with steel reinforcing bars and stirrups using nonlinear 3D finite element (FE) analysis. Test results of 17 large-scale beam specimens tested under monotonic flexural four-point loading from the literature are used as an experimental database to validate the developed nonlinear 3D FE analysis and to study the contributions of steel fibers on the initial stiffness, strength, deformation capacity, cracking behavior, and residual stress. The examined SFRC beams include various ratios of longitudinal reinforcement (0.3%, 0.6%, and 1.0%) and steel fiber volume fractions (from 0.3 to 1.5%). The proposed FE analysis employs the nonlinearities of the materials with new and established constitutive relationships for the SFRC under compression and tension based on experimental data. Especially for the tensional response of SFRC, an efficient smeared crack approach is proposed that utilizes the fracture properties of the material utilizing special stress versus crack width relations with tension softening for the post-cracking SFRC tensile response instead of stress–strain laws. The post-cracking tensile behavior of the SFRC near the reinforcing bars is modeled by a tension stiffening model that considers the SFRC fracture properties, the steel fiber interaction in cracked concrete, and the bond behavior of steel bars. The model validation is carried out comparing the computed key overall and local responses and responses measured in the tests. Extensive comparisons between numerical and experimental results reveal that a reliable and computationally-efficient model captures well the key aspects of the response, such as the SFRC tension softening, the tension stiffening effect, the bending moment–curvature envelope, and the favorable contribution of the steel fibers on the residual response. The results of this study reveal the favorable influence of steel fibers on the flexural behavior, the cracking performance, and the post-cracking residual stress.


Introduction
The addition of randomly distributed discrete fibers in concrete significantly improves the overall performance of reinforced concrete (RC) structural members because of its enhanced tensile properties and cracking control. This improvement is mainly attributed to the crack-bridging phenomenon that is observed across crack surfaces due to the incorporated steel fibers. The advanced characteristics of fiber-reinforced cement-based elements depend on many factors such as size, type, elastic

Constitutive Laws of the Materials and FE Model Formulation
The main aspects of the SFRC analyses developed in the current research are presented in this section. The constitutive laws of SFRC have been derived from experimental results and models proposed by the authors in previous studies available in the literature. Especially for the simulation of the SFRC under tension, a smeared crack approach is adopted. The proposed smeared crack model utilizes the fracture properties of the material and employs constitutive relationships of stress versus crack width with tension softening for the post-cracking SFRC tensile response instead of stress-strain laws. Further, the post-cracking tensile behavior of the SFRC near the reinforcing bars is modeled by a tension stiffening model that considers the SFRC fracture properties and the reinforcement characteristics. Figure 1 illustrates an idealized moment versus curvature response of a SFRC member including four stages of cracking. In the first stage (Stage I), concrete, as well as steel, are considered to behave in a linear and elastic way. Stresses and strains are uniformly distributed along the beam, and the applied load is shared between concrete and reinforcement. Concrete is considered macroscopically uncracked, and the tensile stress in concrete does not exceed its tensile strength. Steel fibers are assumed to not contribute at this stage, as a crack must be formed for the slip between concrete and fiber to occur, and fiber to be activated.

Stages of SFRC Cracking and Tension Stiffening
The second stage (Stage II) begins when the first crack is formed. Once tensile stress overcomes the tensile strength of SFRC and cracking moment is exceeded, a flexural crack begins to develop in the weakest cross-section of the beam. At this stage, the concrete-steel bond mechanism is activated, and the tension stiffening effect starts contributing to the member response. Tension stiffening is the ability of concrete between the cracks to bear tensile stresses, and through the bond mechanism, tensile stresses can be transferred from reinforcement to concrete. In addition, fibers tend to "bridge" the formed cracks, improve the bond efficiency due to the ability to also carry tension across cracks and thereby improve the tension stiffening performance of the member. As a result, the average tensile stiffness of cracked RC is greater than that of the reinforcing bars alone. Therefore, if this effect is disregarded, the stiffness of the cracked RC may be underestimated. As the imposed load increases, more cracks are developed, and the process goes on until the final pattern of cracking is established. If the space between the cracks is not large enough to cultivate a sufficient bond to allow the developed stress to reach tensile strength, no further cracking will occur.
The third stage (Stage III) occurs when the crack pattern has stabilized, and the crack spacing remains constant; cracks self-propagate at this point. The tension stiffening effect decays as the load increases. In the final stage (Stage IV), as the member's average tensile stress begins to exceed the yield

Concrete Damaged Plasticity
The well-known concrete damaged plasticity (CDP) model is usually employed in FE analyses using Abaqus to simulate the nonlinear behavior of concrete. CDP is a general-purpose commonly used model that can be applied to all types of concrete structures under monotonic, cyclic, and dynamic loading. It is a continuum, plasticity-based, damage model that assumes tensile cracking and compressive crushing are the two main failure mechanisms in concrete [37]. In the FE analyses performed in this study, the CDP model has been used for comparison reasons to evaluate the effectiveness and the accuracy of the proposed materials model. Some main aspects and equations of the CDP model are briefly reported below.
The CDP model adopts the elastic model to define the mechanical properties of the materials in the elastic stage. However, when entering the stage of damage, to describe the stiffness degradation, the modulus of elasticity is described as: where E0 is the initial elastic modulus and d is the plastic damage factor (dc or dt, for compression and tension, respectively), which varies from 0 ≤ d ≤ 1, with zero indicating the undamaged material, to one indicating a complete loss of strength. According to Lubliner et al. [37], plastic degradation occurs only within the softening range and the stiffness is proportional to the material's cohesion. Solving the previous equation for the plastic damage factor, it is derived: where c is cohesion in the yield criteria, which is proportional to stress, and cmax is proportional to the strength of the concrete. Adopting Equation (2) to uniaxial tension or compression converts to the following expression: where f is either the compressive or the tensile strength of concrete and σ is the compressive or the tensile stress, respectively.

Concrete Damaged Plasticity
The well-known concrete damaged plasticity (CDP) model is usually employed in FE analyses using Abaqus to simulate the nonlinear behavior of concrete. CDP is a general-purpose commonly used model that can be applied to all types of concrete structures under monotonic, cyclic, and dynamic loading. It is a continuum, plasticity-based, damage model that assumes tensile cracking and compressive crushing are the two main failure mechanisms in concrete [37]. In the FE analyses performed in this study, the CDP model has been used for comparison reasons to evaluate the effectiveness and the accuracy of the proposed materials model. Some main aspects and equations of the CDP model are 5 of 25 where f is either the compressive or the tensile strength of concrete and σ is the compressive or the tensile stress, respectively.
The uniaxial compressive and tensile responses of concrete in relation to the concrete damage plasticity model subjected to compression and tension load are calculated based on the following expressions (see also Figure 2a,b for notation): where σ c − ε c and σ t − ε t are the concrete compressive and tensile stress-strain relationships, respectively, ε c,pl and ε t,pl are the compressive and tensile plastic strain of concrete, respectively.
Materials 2020, 13, x FOR PEER REVIEW 5 of 24 The uniaxial compressive and tensile responses of concrete in relation to the concrete damage plasticity model subjected to compression and tension load are calculated based on the following expressions (see also Figure 2a,b for notation): where − ε and − ε are the concrete compressive and tensile stress-strain relationships, respectively, ε , and ε , are the compressive and tensile plastic strain of concrete, respectively. Using the above equations, the uniaxial stress-strain curves can be converted into stress versus plastic strain curves. This conversion is carried out by Abaqus automatically when user provides in a tabular form all the variables contained in the above equations. The strain variables for the compressive behavior are calculated as follows (see also Figure 2a for notation): The equivalent plastic strain for crushed concrete, εc,pl, is estimated using: where dc is the compressive damage variable: where and σcu = fcm (mean concrete compressive strength).
In the same manner, the cracking strain εt,cr and εt0,el for the tensile behavior are defined as (see also Figure 2b for notation): where the plastic strain for damaged concrete is: with dt the tensile damage variable: Using the above equations, the uniaxial stress-strain curves can be converted into stress versus plastic strain curves. This conversion is carried out by Abaqus automatically when user provides in a tabular form all the variables contained in the above equations. The strain variables for the compressive behavior are calculated as follows (see also Figure 2a for notation): The equivalent plastic strain for crushed concrete, ε c,pl , is estimated using: where d c is the compressive damage variable: where and σ cu = f cm (mean concrete compressive strength).
In the same manner, the cracking strain ε t,cr and ε t0,el for the tensile behavior are defined as (see also Figure 2b for notation): where the plastic strain for damaged concrete is: with d t the tensile damage variable: and σ to = f ctm (mean concrete tensile strength). Five other CDP material associated parameters must be described to completely define the inelastic behavior of concrete [38][39][40]: • ψ is the dilatation angle, which affects the amount of plastic volume deformation. Various values (20 • -45 • ) are found in the literature. A dilation angle close to the friction angle of the material, for concrete 56.3 • , results in a ductile behavior, whereas a low value near 0 • leads to very brittle behavior. • ∈, the flow potential eccentricity, which is the rate of approach of the plastic potential hyperbolic to its asymptote, is used to define the shape of the plastic potential surface in the meridional plane. • σ bo /σ co is the ratio of the strength in the biaxial state to the strength in the uniaxial state. • K c is the ratio of the tensile to the compressive meridian and determines the shape of the yield surface. • µ is the viscosity parameter, which is used to help the analysis achieve a good convergence.
The values used in the current study are presented in Table 1.

Proposed Model of the SFRC Compressive Behavior
Steel fibers practically become effective after cracking and consequently influence post-peak compressive behavior of SFRC. The amount, geometry, and the bond characteristics of the added steel fibers are the parameters affecting the improved post-cracking ductility of SFRC under compression. This enhancement is caused by the ability of steel fibers to exhibit a progressive debonding failure that provides crack-bridging, confinement, and crack growth resisting of the developed cracks [41][42][43][44]. However, compressive strength seems to be slightly increased, especially in SFRC mixtures with a low volume fraction of steel fibers. Most analytical models proposed for the simulation of the SFRC compressive response are based on regression analysis of experimental results available in the literature and taking into account the previously mentioned parameters [45][46][47][48][49][50].
The model adopted in this paper to simulate the compressive stress-strain (σ c -ε c ) behavior of SFRC was proposed by Chalioris and Panagiotopoulos [51], and was derived from test data of 125 stress-strain curves and 257 strength values. Based on this model, the initial ascending compressive behavior until the maximum strength of common SFRC mixtures with cylinder compressive strength, f c,SF , less than 50 MPa is described by the expression (see also Figure 3 for notation): where ε co,SF is the SFRC strain that corresponds to the ultimate compressive stress.
where fc is the cylinder compressive strength of plain concrete, εco is the strain corresponding to the ultimate stress of plain concrete and usually takes the value of 0.002, and F is fiber factor that depends on the characteristics of the added steel fibers in the SFRC mixture and can be evaluated as: = ( ⁄ ), where β is a bond factor that equals to 0.50 for round fibers, 0.75 for deformed (such as hooked, crimped, and undulated) fibers, and 1.0 for indented fibers, and VSF, lSF, and dSF are the volume fraction, length, the diameter of the steel fibers, respectively. The ability of this analytical model to efficiently predict the overall compressive behavior of SFRC was verified against experimental stress-strain curves of various SFRC mixtures by Chalioris and Panagiotopoulos [51].

Proposed Smeared Crack Model of the SFRC Tensile Behavior
Steel fibers significantly improve the tensile performance of SFRC increasing the strength, mainly enhancing the deformation capability of the post-peak behavior by providing pseudo-ductile response due to their gradual debonding procedure [52]. Many analytical stress-strain expressions have been developed to describe the SFRC response under tension [53,54]. In this paper, a smeared crack analysis that has previously been developed and experimentally verified by the authors [55,56] to model the behavior of plain concrete under tension is properly modified and adopted for the SFRC tensile behavior. Similar smeared crack models have also been applied recently in FE simulations to evaluate the uncertainty of crack width of RC beams [57]. The adopted smeared crack model is properly modified herein to employ the favorable contribution of the added steel fibers on the tensile response.
The proposed model employs constitutive relationships in terms of normal stress versus crack width for the post-cracking tensile response of SFRC instead of describing the crack process by stressstrain laws. Based on this aspect, crack propagation of SFRC takes place with the formation of a fracture process zone that is initiated at maximum tensile strength of SFRC, ft,SF, and is characterized by the gradually reduction of the strength during deformation. This fracture process zone is defined as the boundary of the strain softening region, which is considered as a SFRC property and is assumed to be wider than the region of visible cracks. It is also assumed that between the cracks of this zone there are less damaged or even elastic parts. Thus, the total tensional strain, εt, is considered The post-peak response is assumed to be linear from the value of the compressive strength, f c,SF , until the value of 0.85 f c,SF . The parameters of the strength, f c,SF , the strain at ultimate strength, ε co,SF , and the ultimate strain, ε cu,SF , that corresponds to the post-peak compressive stress 0.85 f c,SF are calculated using the following expressions: where f c is the cylinder compressive strength of plain concrete, ε co is the strain corresponding to the ultimate stress of plain concrete and usually takes the value of 0.002, and F is fiber factor that depends on the characteristics of the added steel fibers in the SFRC mixture and can be evaluated as: , where β is a bond factor that equals to 0.50 for round fibers, 0.75 for deformed (such as hooked, crimped, and undulated) fibers, and 1.0 for indented fibers, and V SF , l SF , and d SF are the volume fraction, length, the diameter of the steel fibers, respectively. The ability of this analytical model to efficiently predict the overall compressive behavior of SFRC was verified against experimental stress-strain curves of various SFRC mixtures by Chalioris and Panagiotopoulos [51].

Proposed Smeared Crack Model of the SFRC Tensile Behavior
Steel fibers significantly improve the tensile performance of SFRC increasing the strength, mainly enhancing the deformation capability of the post-peak behavior by providing pseudo-ductile response due to their gradual debonding procedure [52]. Many analytical stress-strain expressions have been developed to describe the SFRC response under tension [53,54]. In this paper, a smeared crack analysis that has previously been developed and experimentally verified by the authors [55,56] to model the behavior of plain concrete under tension is properly modified and adopted for the SFRC tensile behavior. Similar smeared crack models have also been applied recently in FE simulations to evaluate the uncertainty of crack width of RC beams [57]. The adopted smeared crack model is properly modified herein to employ the favorable contribution of the added steel fibers on the tensile response.
The proposed model employs constitutive relationships in terms of normal stress versus crack width for the post-cracking tensile response of SFRC instead of describing the crack process by stress-strain laws. Based on this aspect, crack propagation of SFRC takes place with the formation of a fracture process zone that is initiated at maximum tensile strength of SFRC, f t,SF , and is characterized by the gradually reduction of the strength during deformation. This fracture process zone is defined as the boundary of the strain softening region, which is considered as a SFRC property and is assumed to be wider than the region of visible cracks. It is also assumed that between the cracks of this zone there are less damaged or even elastic parts. Thus, the total tensional strain, ε t , is considered as the sum of an elastic, ε t,el , and a fracture component, ε t,fr , which can be estimated based on the following relationships (see also Figure 4): where σ t is the tensile stress, E t,SF is the modulus of elasticity under tension, w t is the crack width, and L fr,SF is the length of fracture process zone of the SFRC.
Materials 2020, 13, x FOR PEER REVIEW 8 of 24 as the sum of an elastic, εt,el, and a fracture component, εt,fr, which can be estimated based on the following relationships (see also Figure 4): where σt is the tensile stress, Et,SF is the modulus of elasticity under tension, wt is the crack width, and Lfr,SF is the length of fracture process zone of the SFRC.
(a) (b) It is noted that the influence of cracks on the deformations of the beams is considered through the implementation of the aforementioned smeared crack approach to the nonlinear FE analysis using Abaqus software to provide numerical behavioral curves and cracking patterns based on the stress distribution of the simulated SFRC beams.

Elastic Response
The parameters of the elastic component of the proposed model for the simulation of the SFRC tensile behavior, as shown by the σt-εt,el linear diagram in Figure 4, are calculated based on the properties of the corresponding plain concrete and the added steel fibers, as shown in the following relationships [58]: where Et,SF is the modulus of elasticity under tension of SFRC; ft, εto, and Et are the ultimate tensile strength, the corresponding strain, and the modulus of elasticity under tension of the plain concrete, respectively; fSF, εy,SF, and ESF are the ultimate tensile strength at yielding, the yield strain, and the modulus of elasticity of the steel fiber, respectively; and nl,el is the ratio of the average fiber elastic stress to the maximum fiber stress that is usually equal to 0.5. It is noted that the influence of cracks on the deformations of the beams is considered through the implementation of the aforementioned smeared crack approach to the nonlinear FE analysis using Abaqus software to provide numerical behavioral curves and cracking patterns based on the stress distribution of the simulated SFRC beams.

Elastic Response
The parameters of the elastic component of the proposed model for the simulation of the SFRC tensile behavior, as shown by the σ t -ε t,el linear diagram in Figure 4, are calculated based on the properties of the corresponding plain concrete and the added steel fibers, as shown in the following relationships [58]: where E t,SF is the modulus of elasticity under tension of SFRC; f t , ε to , and E t are the ultimate tensile strength, the corresponding strain, and the modulus of elasticity under tension of the plain concrete, respectively; f SF , ε y,SF , and E SF are the ultimate tensile strength at yielding, the yield strain, and the modulus of elasticity of the steel fiber, respectively; and n l,el is the ratio of the average fiber elastic stress to the maximum fiber stress that is usually equal to 0.5.

Fracture Response with Tension Softening
The fracture characteristics and the tension softening response of SFRC define the parameters of the fracture component of the proposed smeared crack model. The energy required for the formation of the cracks included in the fracture process zone and for the fully opening of one single crack for a unit area crack plane is the fracture energy, G f,SF , and can be expressed as: The post-peak fracture response (σ t -w t curve in Figure 4) is described by a linear descending line until the point of the maximum post-cracking stress, k f f t,SF , and the corresponding crack width, k w w u,SF , and then a horizontal line of constant tensile stress (σ t = k t f t,SF when w t > k w w u,SF ) until the ultimate considered crack width, w u,SF . Thus, the fracture energy can also be expressed in terms of the area under the curve of SFRC tensile stress versus crack width, as shown by the σ t -w t bilinear diagram in Figure 4: The values of the maximum post-cracking stress and crack width depend on the SFRC characteristics and can be estimated using the following expressions for short deformed steel fibers [58]: where σ fu is the ultimate stress of the fiber when a uniform ultimate bond stress, τ u , is assumed at the fiber-matrix interface and calculated as follows: where l cr is the length in the fiber required to develop the ultimate fiber stress that can be estimated as:

Fracture Energy
The value of the SFRC fracture energy is an important parameter that depends on the SFRC characteristics and can be evaluated experimentally by tension tests. In this study, the ratio of the fracture energy of the SFRC, G f,SF , to the fracture energy of the corresponding plain concrete, G f , was determined by 85 test results from the literature as a function of the fiber factor, F, as shown in Figure 5. The experimental results of low to high concrete strength deformed (hooked, crimped, and undulated) short steel fibers were used [59][60][61][62][63][64]. The experimental data points and fitted lines in Figure 5 include test results of low concrete strength with steel fibers of normal strength (Line A), normal to high concrete strength with steel fibers of high strength (Line B), and normal to high concrete strength with steel fibers of normal strength (Line C). Thus, the expression derived from the experimental results presented in Figure 5 has the form: where kG is a factor that can be taken based on the results in Figure 5. The fracture energy of the plain concrete, Gf, can be expressed using the following relationship assuming a linear reduction of the tensile concrete stress from the maximum tensile strength, ft, to zero at the ultimate crack width, wu: where afr is a coefficient that depends on the shape of the stress versus crack width curve and the nature and size of the concrete aggregates and takes values from 5 to 8 for concrete aggregates of maximum size dg = 32 to 8 mm, respectively [55]. Lfr is the length of fracture process zone of plain concrete that can be taken as 3dg [65]. Thus, using Equations (33)- (35), Equation (27) can be written as: , , Based on the aspects of the proposed smeared crack model, the formulation of the stress is:

Fracture Process Zone of the SFRC
The aforementioned formulation can provide a feasible evaluation of the entire tensile behavior of the SFRC based on the elastic and fracture characteristics of the materials used in the fibrous mixture (concrete and steel fibers). All parameters and coefficients can be defined using the tensile strength, ft, the modulus of elasticity under tension, Et, and the maximum size of the aggregates, dg, of the plain concrete, along with the characteristics of the added steel fibers, such as their ultimate tensile strength at yielding, fSF, their modulus of elasticity, ESF, and the fiber factor, F, which depend on their characteristics (bond factor, β, volume fraction, VSF, length, lSF, and diameter, dSF). The only  Thus, the expression derived from the experimental results presented in Figure 5 has the form: where k G is a factor that can be taken based on the results in Figure 5. The fracture energy of the plain concrete, G f , can be expressed using the following relationship assuming a linear reduction of the tensile concrete stress from the maximum tensile strength, f t , to zero at the ultimate crack width, w u : where a fr is a coefficient that depends on the shape of the stress versus crack width curve and the nature and size of the concrete aggregates and takes values from 5 to 8 for concrete aggregates of maximum size d g = 32 to 8 mm, respectively [55]. L fr is the length of fracture process zone of plain concrete that can be taken as 3d g [65]. Thus, using Equations (33)- (35), Equation (27) can be written as: Based on the aspects of the proposed smeared crack model, the formulation of the stress is:

Fracture Process Zone of the SFRC
The aforementioned formulation can provide a feasible evaluation of the entire tensile behavior of the SFRC based on the elastic and fracture characteristics of the materials used in the fibrous mixture (concrete and steel fibers). All parameters and coefficients can be defined using the tensile strength, f t , the modulus of elasticity under tension, E t , and the maximum size of the aggregates, d g , of the plain concrete, along with the characteristics of the added steel fibers, such as their ultimate tensile strength at yielding, f SF , their modulus of elasticity, E SF , and the fiber factor, F, which depend on their characteristics (bond factor, β, volume fraction, V SF , length, l SF , and diameter, d SF ). The only parameter that requires being defined is the size of the SFRC fracture process zone length, L fr,SF . This length cannot be assumed as mesh dependent, but it has to be attributed to the nature of the SFRC, the type of loading and the size of the specimen [55].
In this study the estimation of the crack process zone length is based on experimental data and the developed FE modeling. In the following sections, the methodology used to estimate the value of the fracture process zone length of the SFRC is presented.

Modeling of Steel Reinforcement
Steel reinforcement (longitudinal and stirrups) was modeled as a linear elastic material until the reach of yield point. Post-yielding behavior is considered to be nonlinear inelastic, referred as plastic behavior. Plastic behavior is characterized by permanent (plastic) deformations and is defined by introducing parameters of stress and plastic strain at the yield point (f y , ε y ) as well as the stress and strain at ultimate point (f y and ε u ). Steel modulus of elasticity, E s , is also provided to Abaqus, according to experimental data values as well as Poisson's ratio, which is considered equal to 0.3 for all analyses performed in this study.

Examined Beam Specimens
The current research presents a tension stiffening FE model for SFRC that properly considers the above-mentioned parameters. Seventeen SFRC beams available in the literature [66][67][68][69][70][71] were selected and a nonlinear FE analysis was performed to illustrate the capability and accuracy of the proposed model. The selected series of experimental tests had the same geometrical characteristics but differentiated in longitudinal reinforcement, concrete strength, and amount of steel fibers added to the concrete mixture.
These test series were chosen because all the necessary test data needed by the FE simulations were clearly reported. Further, the beam specimens include various ratios of longitudinal reinforcement (0.3%, 0.6%, and 1.0%) and steel fiber volume fractions (from 0.3 to 1.5%). Table 2 summarizes the geometric dimensions and reinforcement details of all examined beams as well as the material properties.
The 17 FE simulations were developed using same geometry, dimensions, material properties, and boundary conditions of the tested beams. The model is based on a 3D FE analysis implemented in Abaqus [34] software and on a smeared cracked approach, which does not require predetermined crack paths. Proper modeling and analysis of the SFRC beams involves full consideration on the material and geometry nonlinearities as well as effects related to the interaction and load sharing between steel and reinforced concrete. More details about the FE model are presented in the following sections.

Element Types Selected
SFRC was modeled using three-dimensional eight-node solid elements with reduced integration to prevent shear locking effect (C3D8R). For steel reinforcement (bars and stirrups), three-dimensional two-node truss elements were chosen (T3D2). Each node of these elements has three degrees of freedom with translation in X, Y, and Z directions (global coordinates system), as shown in Figure 6.

Element Types Selected
SFRC was modeled using three-dimensional eight-node solid elements with reduced integration to prevent shear locking effect (C3D8R). For steel reinforcement (bars and stirrups), threedimensional two-node truss elements were chosen (T3D2). Each node of these elements has three degrees of freedom with translation in X, Y, and Z directions (global coordinates system), as shown in Figure 6. Proper modeling of steel reinforcement requires that bond interaction between reinforcement and concrete is considered. In general, the steel reinforcement could be modeled using a beam element with bending stiffness. However, in this study, the bending stiffness of the reinforcement is considered to be quite low compared to the bending stiffness of the SFRC matrix, and therefore it is neglected. Thus, truss elements are used for modeling steel reinforcement. The interaction between concrete and reinforcement is considered using the embedded method and specifically "the truss in solid" type. Truss elements (longitudinal reinforcement and stirrups) are embedded in solid region (SFRC), referred as host region. Geometric relationships between host and embedded elements are explored and if an embedded element node lies within a host element, then the transactional degrees of freedom (DOFs) of the node are constrained to follow the response of the host element.

Boundary Conditions and Loading
Displacement boundary conditions are necessary to constrain the model and obtain a unique solution. To ensure that the model behaves in an equivalent way as the experimental beam, the supports and loadings must be introduced to proper boundary conditions. The entire SFRC beam was modeled and simply supported. Two supports were designed in such a way as to create one roller and one pinned support. The supports were placed at a distance of 140 mm from each end of the beam according to the experimental set up, while the ends of the beams were free. At the one support, in the Ux, Uy, and Uz directions, a single line of nodes was constrained and added as constant values of 0, while at the other support only Uy direction were restricted given value of 0, as shown in Figure 7. This allows the beam to rotate at the supports. Proper modeling of steel reinforcement requires that bond interaction between reinforcement and concrete is considered. In general, the steel reinforcement could be modeled using a beam element with bending stiffness. However, in this study, the bending stiffness of the reinforcement is considered to be quite low compared to the bending stiffness of the SFRC matrix, and therefore it is neglected. Thus, truss elements are used for modeling steel reinforcement. The interaction between concrete and reinforcement is considered using the embedded method and specifically "the truss in solid" type. Truss elements (longitudinal reinforcement and stirrups) are embedded in solid region (SFRC), referred as host region. Geometric relationships between host and embedded elements are explored and if an embedded element node lies within a host element, then the transactional degrees of freedom (DOFs) of the node are constrained to follow the response of the host element.

Boundary Conditions and Loading
Displacement boundary conditions are necessary to constrain the model and obtain a unique solution. To ensure that the model behaves in an equivalent way as the experimental beam, the supports and loadings must be introduced to proper boundary conditions. The entire SFRC beam was modeled and simply supported. Two supports were designed in such a way as to create one roller and one pinned support. The supports were placed at a distance of 140 mm from each end of the beam according to the experimental set up, while the ends of the beams were free. At the one support, in the Ux, Uy, and Uz directions, a single line of nodes was constrained and added as constant values of 0, while at the other support only Uy direction were restricted given value of 0, as shown in Figure 7. This allows the beam to rotate at the supports. The load was applied as a distributed surface load (Figure 7), which is more realistic as in real structures forces are most likely to be transferred through contact and thus a force is likely to be distributed across a contact area. Applying the load as a concentrated force to a single node can cause troublesome effects, such as high stress concentration around the application node. When a load is applied as surface load, the total force is distributed to the face of the FE which composes the surface, The load was applied as a distributed surface load (Figure 7), which is more realistic as in real structures forces are most likely to be transferred through contact and thus a force is likely to be distributed across a contact area. Applying the load as a concentrated force to a single node can cause troublesome effects, such as high stress concentration around the application node. When a load is applied as surface load, the total force is distributed to the face of the FE which composes the surface, and then to the nodes laying on the face, which leads to having a smaller amount of the total force at each node and excessive local stress concentration is avoided. Further, to ensure that the dynamic implicit procedure adopted in this study obtained a quasi-static solution, the load was applied steadily and gradually to avoid any significant change in acceleration from one increment to the next. This further guarantees that the alterations in stresses and displacement are smooth. The loading sequence used in the FE analysis was based on the experimental sequence and the analysis was performed using the load control technique.

Meshing and Convergence
The mesh size of 50 mm was applied to all elements (truss and solid) to make sure that different materials (steel and concrete) share common nodes. This size of mesh was chosen because concrete cracking propagation usually involves spatial scales in the range of 2-3 dominant aggregate sizes (usually between 16 and 32 mm) of the concrete [65]. A convergence study was carried out using different mesh sizes before proceeding to the main analysis and the above-mentioned mesh proved to be fine enough to generate reliable results.

Material Input
The main properties of the proposed material models for each beam are presented in Table 2. Some additional parameters also used as input in the performed nonlinear analysis of the SFRC beams are presented in Table 3. Stress and crack widths at Points 1-3 and corresponding plastic damage factors shown in Table 3 are defined in the SFRC tension model depicted in Figure 4. Table 3. Input parameters of the examined SFRC beams.

Ref.
Beam Name

Validation and Accuracy of the Proposed Model
As mentioned in Section 2.4.4, the value of the SFRC fracture process zone length, L fr,SF , cannot be assumed as mesh dependent but it has to be attributed to the nature of the SFRC, the type of loading, and the size of the specimen. In this study the estimation of the proper length size for the specific loading case is based on the examined tests using the developed FE modeling and by the comparison of the numerical results and the experimental data.
Thus, comparisons of the numerically predicted bending moment versus curvature responses of SFRC beams S3-1-F05 shown in Figure 8 for different values of fracture process zone length, L fr,SF , with the experimental response can yield useful conclusions about the proper size of this length and for this type of loading. The nature of the material used in this procedure is represented by the adopted stress versus crack width curve. The examined values of the fracture process zone length, L fr,SF , are taken as function of the steel fiber length, l SF , and five different values are investigated: L fr,SF = 2.0, 2.5, 3.0, 3.5 and 4.0 l SF . The comparisons of these five numerical curves (blue continuous lines) with the experimental data (red dotted line) shown in Figure 8 indicate that the best fitted numerical curve has been derived for L fr,SF = 3l SF . Further, the comparisons between the numerical predictions using this specific length of the SFRC fracture process zone the test results prove that the proposed model fits very well to the experimental behavior of the beams with higher accuracy than the simulation without considering the tension stiffening and softening effect.
As mentioned in Section 2.4.4, the value of the SFRC fracture process zone length, Lfr,SF, cannot be assumed as mesh dependent but it has to be attributed to the nature of the SFRC, the type of loading, and the size of the specimen. In this study the estimation of the proper length size for the specific loading case is based on the examined tests using the developed FE modeling and by the comparison of the numerical results and the experimental data.
Thus, comparisons of the numerically predicted bending moment versus curvature responses of SFRC beams S3-1-F05 shown in Figure 8 for different values of fracture process zone length, Lfr,SF, with the experimental response can yield useful conclusions about the proper size of this length and for this type of loading. The nature of the material used in this procedure is represented by the adopted stress versus crack width curve. The examined values of the fracture process zone length, Lfr,SF, are taken as function of the steel fiber length, lSF, and five different values are investigated: Lfr,SF = 2.0, 2.5, 3.0, 3.5 and 4.0 lSF. The comparisons of these five numerical curves (blue continuous lines) with the experimental data (red dotted line) shown in Figure 8 indicate that the best fitted numerical curve has been derived for Lfr,SF = 3lSF. Further, the comparisons between the numerical predictions using this specific length of the SFRC fracture process zone the test results prove that the proposed model fits very well to the experimental behavior of the beams with higher accuracy than the simulation without considering the tension stiffening and softening effect.   [67,68] with 0.5%, 1.0%, 1.5%, and 1.0% volume fraction of steel fibers, respectively. To comprehend the influence of the proposed smeared crack model with tension softening and the tension stiffening effect on the flexural behavior of the examined beams, two curves derived from FE analysis, in terms of bending moment versus curvature are compared with the experimental one in Figure 9a, Figure 10a, Figure 11a, and Figure 12a. The first curve (blue continuous line) was derived using the proposed material models in the nonlinear analysis with tension softening and tension stiffening approach and the second (black thin dotted line) using the FE model without considering these effects (denoted as "model without TS" where TS is tension softening and stiffening).
The presented stress versus strain curves in Figures 9b, 10b, 11b and 12b consist of the stresses corresponding to the tension stiffening effect and the residual stresses due to fiber interaction with the concrete [67,68]. The numerical curves in these diagrams were calculated based on the stress-strain values in each element, which were extracted by the FE analysis and then the average values over all elements were calculated and used. The averaged stress strain values are the overall stress versus strain for the whole beam simulation and combined in a stress-strain curve.
behavior of the examined beams, two curves derived from FE analysis, in terms of bending moment versus curvature are compared with the experimental one in Figures 9a, 10a, 11a, and 12a. The first curve (blue continuous line) was derived using the proposed material models in the nonlinear analysis with tension softening and tension stiffening approach and the second (black thin dotted line) using the FE model without considering these effects (denoted as "model without TS" where TS is tension softening and stiffening). proposed smeared crack model with tension softening and the tension stiffening effect on the flexural behavior of the examined beams, two curves derived from FE analysis, in terms of bending moment versus curvature are compared with the experimental one in Figures 9a, 10a, 11a, and 12a. The first curve (blue continuous line) was derived using the proposed material models in the nonlinear analysis with tension softening and tension stiffening approach and the second (black thin dotted line) using the FE model without considering these effects (denoted as "model without TS" where TS is tension softening and stiffening).    The presented stress versus strain curves in Figures 9b, 10b, 11b and 12b consist of the stresses corresponding to the tension stiffening effect and the residual stresses due to fiber interaction with the concrete [67,68]. The numerical curves in these diagrams were calculated based on the stressstrain values in each element, which were extracted by the FE analysis and then the average values over all elements were calculated and used. The averaged stress strain values are the overall stress versus strain for the whole beam simulation and combined in a stress-strain curve.
It is noted that the point of the peak stress in the stress-strain diagrams of Figures 9b, 10b, 11b, and 12b corresponds to the cracking point in the moment versus curvature diagram shown in Figures  9a, 10a, 11a, and 12a, respectively. Since steel fibers mainly influence the post-cracking behavior of SFRC, the post-peak descending part of the residual stress versus strain curve depends on the volume fraction of the added steel fibers. Thus, SFRC mixtures with higher amount of fibers (for example, beam S3-1-F15 with VSF = 1.5%) exhibit higher post-cracking stress.
Further, the numerical predictions of the stress versus strain curves using the proposed nonlinear smeared crack model with tension softening and tension stiffening are in very good Numerical curves derived by the developed nonlinear FE model with and without tension softening and stiffening effect are also compared with the experimental ones in the diagrams of Figure 13 for the SFRC beams: (a) S3-2-7F; (b) B2-F03; (c) B3-F06; (d) S1-F05; (e) S1-F10; and (f) S1-F15 [66,70,71]. These comparisons reveal that in most of the examined cases the numerical predictions of the proposed FE model that considers the smeared crack analysis with tension softening and tension stiffening effect are in closer agreement with the test data than the numerical curves without this effect. Numerical curves derived by the developed nonlinear FE model with and without tension softening and stiffening effect are also compared with the experimental ones in the diagrams of Figure 13 for the SFRC beams: (a) S3-2-7F; (b) B2-F03; (c) B3-F06; (d) S1-F05; (e) S1-F10; and (f) S1-F15 [66,70,71]. These comparisons reveal that in most of the examined cases the numerical predictions of the proposed FE model that considers the smeared crack analysis with tension softening and tension stiffening effect are in closer agreement with the test data than the numerical curves without this effect. The effectiveness of the proposed nonlinear FE model to predict accurately the flexural response of SFRC beams with different longitudinal reinforcement ratios is examined in the diagrams of Figure  14 in terms of bending moment versus curvature behavioral curves. The SFRC beams compared in the diagrams of Figure 14a-c has steel fiber volume fraction 0.5%, 1.0%, and 1.5%, respectively. Further, the ratio of the steel longitudinal reinforcement of the beams are as follows: • Beams of series "S1" (S1-F05, S1-F10 and S1-F15): ρl = 1.0%. The effectiveness of the proposed nonlinear FE model to predict accurately the flexural response of SFRC beams with different longitudinal reinforcement ratios is examined in the diagrams of Figure 14 in terms of bending moment versus curvature behavioral curves. The SFRC beams compared in the diagrams of Figure 14a-c has steel fiber volume fraction 0.5%, 1.0%, and 1.5%, respectively. Further, the ratio of the steel longitudinal reinforcement of the beams are as follows: • Beams of series "S1" (S1-F05, S1-F10 and S1-F15): ρ l = 1.0%. The comparisons shown in Figure 14 reveal that a good agreement between the numerical and the experimental bending moment versus curvature curves was achieved. The comparisons shown in Figure 14 reveal that a good agreement between the numerical and the experimental bending moment versus curvature curves was achieved. Further, to establish the validity of the proposed nonlinear FE analysis, Table 4 presents the calculated values of the differences between the numerical predictions and the test results in terms of "calculation errors". Specifically, the discrepancy of the bending moment values between the numerical predictions and the tests for the same curvatures along the entire bending moment versus curvature diagrams of each examined beam was calculated to validate, statistically compare, and Further, to establish the validity of the proposed nonlinear FE analysis, Table 4 presents the calculated values of the differences between the numerical predictions and the test results in terms of "calculation errors". Specifically, the discrepancy of the bending moment values between the numerical predictions and the tests for the same curvatures along the entire bending moment versus curvature diagrams of each examined beam was calculated to validate, statistically compare, and estimate the overall accuracy of the proposed model. This discrepancy in terms of error percentage for each point was calculated according to the following known expression: where M FE model and M exp are the bending moment values derived from the proposed FE model and the experiments, respectively. Table 4 presents the values of the mean absolute error (MAE), the standard error (SE), and coefficient of variation (CV). These values indicate that model predictions lead to an error below 9% for every SFRC beam examined and the average value of MAE of all beams is 5.1%. Figure 15 includes a diagram with the experimental and the numerical behavioral curves of beams S2-F05, S2-F10, and S2-F15 [69]. In the same figure, the schemes of the cracking pattern and the stress distribution of the beams derived by the proposed model at bending moment 50 kNm are also compared. For simplicity and comparison reasons, these numerical schemes show the plane y-z view of the 3D FE model of the beams and only beam S2-F15 is illustrated in both 3D and plane views. It is noted that SFRC beam S2-F15 with higher amount of steel fibers (V SF = 1.5%) demonstrates higher strength and lower deformation than the other SFRC beams with 1.0% and 0.5% steel fiber volume fraction. Further, this SFRC beam (specimen S2-F15) exhibits increased number of cracks and decreased crack spacing and width with respect to the cracking performance of beam S2-F05 with lower amount of steel fibers (V SF = 0.5%) due to the favorable contribution of the steel fibers. These observations are based on the numerical cracking pattern and stress distribution derived from the proposed FE analysis ( Figure 15). The same observations have also been established by relevant tests of flexural SFRC beams [72][73][74].

Influence of Steel Fibers on the Flexural Performance
fraction. Further, this SFRC beam (specimen S2-F15) exhibits increased number of cracks and decreased crack spacing and width with respect to the cracking performance of beam S2-F05 with lower amount of steel fibers (VSF = 0.5%) due to the favorable contribution of the steel fibers. These observations are based on the numerical cracking pattern and stress distribution derived from the proposed FE analysis (Figure 15). The same observations have also been established by relevant tests of flexural SFRC beams [72][73][74]. Further, the contribution of steel fibers on the increase of the residual strength and the overall post-cracking stress versus strain behavior of SFRC structural members with steel bars and stirrups is demonstrated in Figure 16. Figure 16a presents the residual stress versus strain response derived from the proposed model of beams S3-1-F05, S3-1-F10, and S3-1-F15 with VSF = 0.5%, 1.0%, and 1.5% (beams of series "S3-1") and Figure 16b presents the corresponding numerical response of beams S2-F05, S2-F10, and S2-F15 with VSF = 0.5%, 1.0%, and 1.5% (beams of series "S2"). In both examined cases, it is indicated that the increase of the volume fraction of the added steel fibers causes a gradual increase of the residual strength. Further, the contribution of steel fibers on the increase of the residual strength and the overall post-cracking stress versus strain behavior of SFRC structural members with steel bars and stirrups is demonstrated in Figure 16. Figure 16a presents the residual stress versus strain response derived from the proposed model of beams S3-1-F05, S3-1-F10, and S3-1-F15 with V SF = 0.5%, 1.0%, and 1.5% (beams of series "S3-1") and Figure 16b presents the corresponding numerical response of beams S2-F05, S2-F10, and S2-F15 with V SF = 0.5%, 1.0%, and 1.5% (beams of series "S2"). In both examined cases, it is indicated that the increase of the volume fraction of the added steel fibers causes a gradual increase of the residual strength.

Conclusions
The influence of steel fibers on the short-term flexural behavior of realistic SFRC structural members reinforced with steel reinforcement was investigated by nonlinear 3D finite element (FE) analysis using software Abaqus. A database of 17 large-scale beam specimens collected from the literature was used as experimental baseline to perform the developed analysis that considers the tension stiffening effect and the nonlinearities of the materials by new and established constitutive relationships for steel fiber-reinforced concrete (SFRC) under compression and tension. The

Conclusions
The influence of steel fibers on the short-term flexural behavior of realistic SFRC structural members reinforced with steel reinforcement was investigated by nonlinear 3D finite element (FE) analysis using software Abaqus. A database of 17 large-scale beam specimens collected from the literature was used as experimental baseline to perform the developed analysis that considers the tension stiffening effect and the nonlinearities of the materials by new and established constitutive relationships for steel fiber-reinforced concrete (SFRC) under compression and tension. The appropriate loading and constraint conditions were also considered in order for the performed FE analyses to be very close to experimental conditions. Based on the results of this study, the following conclusions can be drawn: • The proposed FE model considers the contribution of the steel fibers according to their type, aspect ratio, and volume fraction to achieve a more realistic evaluation of the compressive and tensile behavior of SFRC utilizing well-established constitutive laws based on experimental data. Especially, for the simulation of the SFRC under tension a smeared crack approach is proposed that utilizes the fracture properties of the material and employs constitutive relationships of stress versus crack width with tension softening for the post-cracking SFRC tensile response instead of stress-strain laws. The post-cracking tensile behavior of the SFRC near the reinforcing bars is modeled by a tension stiffening model that considers the SFRC fracture properties and the reinforcement characteristics.

•
The size of the fracture process zone of the SFRC under tension is attributed to the material, the type of loading, and the size of the specimen. An estimation of this length based on comparisons between numerical and experimental data is presented.

•
The examined realistic SFRC structural members with increased values of the fiber factor, F, exhibit enhanced flexural performance in terms of initial stiffness, strength, deformation capacity, cracking behavior, and residual stress. The developed FE analysis predicts accurately the overall experimental flexural behavior and points out the contribution of the steel fibers on this improvement.

•
Beams with higher amount of steel fibers (V SF = 1.5%) demonstrated lower deformation than the other SFRC beams with 1.0% and 0.5% steel fiber volume fraction at the same level of the applied load. Further, the favorable contribution of the added steel fibers on the cracking performance of the examined SFRC beams is also highlighted in the numerical results of this study. Beams with 1.5% steel fiber volume fraction showed increased number of cracks and decreased crack spacing and width with respect to the corresponding beams with lower amount of fibers (V SF = 0.5%).

•
The favorable contribution of the added steel fibers on the post-cracking short-term behavior of realistic SFRC structural members with reinforcing bars and stirrups was revealed by the performed analyses. Numerical results of the examined SFRC beams show that the post-peak descending part of the residual stress versus strain curves depends on the volume fraction of the added steel fibers and SFRC mixtures with higher amounts of fibers exhibit higher post-cracking stress.

•
Comparisons between experimental results and numerical predictions include typical bending moment versus curvature behavioral curves, deformation-cracking patterns, and residual stress versus strain diagrams. The proposed nonlinear smeared crack analysis with tension softening and tension stiffening effect was found to be reliable and capable of accurately simulating the flexural response of large-scale SFRC beams with steel reinforcement since the predicted behavioral curves are in good agreement with the corresponding experimental ones.
Author Contributions: All authors contributed extensively to this study, discussed the results and reviews, prepared the manuscript, and agreed to the amendments at all stages of the paper. V.K.K. developed the aspects of the numerical investigation and performed the FE analyses under the supervision of C.E.C. and C.G.K. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.