A Numerical Model for Tracing Structural Response of Ultra-High Performance Concrete Beams

This paper presents a finite element-based numerical model for tracing the behavior of ultra-high performance concrete (UHPC) beams. The model developed in ABAQUS can account for stress–strain response of UHPC and reinforcing bar in both tension and compression, bond between concrete and reinforcing steel, and strain hardening effects in bars and UHPC and can trace the detailed response of UHPC beams in the entire range of loading. This model is validated by comparing predicted response parameters including load-strain, load-deflection, and crack propagation against experimental data governed from tests on UHPC beams with different reinforcement ratios, fiber volume fractions, and loading configurations (shear and flexural loading). The validated model is applied to quantify the contribution of stirrups and concrete to shear strength of beams so as to explore the feasibility of removing shear reinforcement in UHPC beams.


Introduction
Ultra-high performance concrete (UHPC) is a new class of cementitious material possessing excellent mechanical properties [1,2]. One of the major drawbacks of UHPC is the brittle behavior of such concretes. This drawback is overcome to certain extent through adding steel fibers to UHPC mix and this type of concrete is referred to as ultra-high performance fiber-reinforced concrete (UHPFRC). UHPFRC exhibits superior strength, improved fracture toughness, ductility, energy absorption capacity, and enhanced post cracking tensile response [2,3]. These outstanding properties result from extremely homogenized microstructure achieved through optimizing the granular mixture along with a low water-to binder ratio, high fineness admixtures, and effective use of steel fibers [4].
In recent years, UHPFRC is finding increasing applications in infrastructural systems owing to its superior properties [5][6][7]. Despite recent research efforts to evaluate structural behavior of UHPC (or UHPFRC), currently there are limited design provisions for structural design of UHPC in international codes and standards (AFGC-SETRA [8], JSCE [9], and KCI [10]).
A number of experiments have been conducted in literature on response of UHPFRC beams. These studies revealed that increasing steel fiber fraction with high specific surface area, improves post cracking stiffness and thus enhances flexural capacity of UHPFRC beams. The results also demonstrate that using a higher fiber volume fraction in UHPFRC and a lower shear span to depth ratio can result in a higher shear capacity [5,[11][12][13][14][15][16][17].
A number of detailed numerical studies has been reported in literature on the response of beams made with different types of concrete such as normal strength concrete (NSC), high strength concrete (HSC), and fiber-reinforced concrete (FRC) [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. There are a number of numerical studies on the structural behavior of UHPFRC beams. Majority of the conducted numerical studies focused on sectional analysis approach to trace moment curvature behavior of beams subjected to flexural loading [5,11]. There are limited finite Modelling 2021, 2 449 element-based numerical studies that evaluated the behavior of UHPFRC members. Tysmans et al. [33] evaluated the behavior of high performance concrete under biaxial tension. These authors inferred that the developed model incorporating concrete damage plasticity (CDP) material model can account for strain hardening behavior of UHPFRC in tension and thus can predict realistic load capacity of UHPFRC member. However, in their model validation, the authors relied upon small-scale specimens with no reinforcing bars. Singh et al. [34] developed a finite element model for evaluating structural behavior of UHPFRC beams wherein developed model was validated only under flexural and not under dominant shear loading. Chen and Graybeal [35,36] concentrated on evaluating load-deflection and strain response of UHPFRC I-girders and Pi-girders subjected to shear and flexure. The authors demonstrated that finite element model applying smeared cracking model produce a stiffer response as compared to CDP model. However, in this study strain hardening response of UHPFRC in tension was neglected. Bahij et al. [37] developed a numerical model for tracing shear response of UHPFRC beams. However, strain hardening in steel and UHPFRC was not included in the developed model and therefore post-peak response of the beams was not captured. Moreover, shear contribution of concrete and stirrups was not quantified.
The developed numerical models were not validated for tracing shear behavior of UHPFRC beams. Further, majority of previous studies focused on global response of UHPFRC structural members with no attention to local response (such as crack propagation direction, contribution of concrete and stirrups to shear capacity) and strain hardening in UHPFRC under tension was neglected. To address some of the noted issues, Solhmirzaei and Kodur [38] developed a finite element numerical model for tracing structural response of UHPFRC beams. However, the model was developed through load control technique and was not able to trace post peak response of the beams. The model was not validated fully and shear response of beams, failure mode (crack propagation direction) were not studied in detail.
This study presents details on the development of a numerical model in tracing the comprehensive structural response of UHPC and UHPFRC beams with a focus on material models to be adopted for evaluating realistic response in the entire range of loading till collapse of member. The model is applied to quantify contribution of concrete and steel to shear strength of UHPFRC beams and to study feasibility of removing shear reinforcement in flexural UHPFRC members.

Finite Element Model
A numerical model for tracing the structural behavior of a UHPFRC beam in the entire range of loading, from preloading stage till collapse, is developed in ABAQUS. A displacement controlled technique is utilized to trace the post peak response of the beam wherein displacement incremented at the nodes located under load points in steps till failure is attained. Details of the developed numerical model are described in following sections.

Discretization of the Beam
In a given UHPC (or UHPFRC) beam, concrete mass is discretized with brick elements (C3D8) while reinforcing steel is modeled as link elements (T3D2). C3D8 element is of eight nodes with three degrees of freedom in directions of x, y, and z. The steel reinforcing bars were modeled using truss elements (T3D2), with the cross section of each bar defined within 3D solid. T3D2 are two noded elements and are used to model one-dimensional reinforcing bars that are assumed to deform by axial stretching only. The bond between concrete and reinforcement is modeled by using bond-link element approach [39,40]. Figure 1 illustrates discretization of a typical UHPFRC beam. Details on size of the beams and reinforcement are presented in Section 3.1.
ing 2021, 2, FOR PEER REVIEW proach [39,40]. Figure 1 illustrates discretization of a typical UHPFRC beam. Detai size of the beams and reinforcement are presented in Section 3.1.

Material Models for UHPC
A damage-based concrete plasticity model (CDP) in ABAQUS was used to trac nonlinear material response of UHPC (or UHPFRC). CDP is a smeared crack ma model and is based on the theory of plastic flow [33]. The yield surface in CDP mo based on the yield surface proposed by Lubliner et al. [41] along with modifica proposed by Lee and Fenves [42] to account for different evolution laws of the stre under tension and compression. In order to define CDP model, a set of material pr ties such as tension stiffening in concrete, response under compression, elastic mod density, Poisson's ratio, and plasticity parameters requires to be input in m [33,34,39,43].

Parameters for Modeling UHPC Behavior
Five parameters relating to plasticity response in UHPC, namely 0 0 ⁄ , , are required to define CDP model. Two parameters of 0 0 ⁄ and modify the surface. 0 0 ⁄ is the ratio of biaxial compressive strength to uniaxial compre strength which influences the yield surface in a plane stress state and the parameter the ratio between distances measured from the hydrostatic axis to tensile and com sive meridians and is used to define the shape of the failure surface in deviatoric p The other two parameters, and modify the non-associated potential flow. is th lation angle that describes the angle of inclination of the failure surface toward th drostatic axis measured in the meridional plane. is an eccentricity parameter that trols the deviation of the hyperbolic plastic potential from its asymptote. is the vi ty parameter which is used for the visco-plastic regularization of the concrete con tive equations [33,35,39,43]. These five parameters were selected to be 1.05, 2/3, 30 1E-4, respectively, according to the literature and sensitivity analysis conducted in study [33,35,38,[43][44][45]. In addition, Poisson's ratio of UHPC (or UHPFRC) were co ered to be 0.2 according to AFGC [8]. Density of UHPC (or UHPFRC) was measur be 2565 kg m 3 ⁄ [17]. CDP model accounts for stiffness degradation in compression and tension.

Material Models for UHPC
A damage-based concrete plasticity model (CDP) in ABAQUS was used to trace the nonlinear material response of UHPC (or UHPFRC). CDP is a smeared crack material model and is based on the theory of plastic flow [33]. The yield surface in CDP model is based on the yield surface proposed by Lubliner et al. [41] along with modifications proposed by Lee and Fenves [42] to account for different evolution laws of the strength under tension and compression. In order to define CDP model, a set of material properties such as tension stiffening in concrete, response under compression, elastic modulus, density, Poisson's ratio, and plasticity parameters requires to be input in model [33,34,39,43].

Parameters for Modeling UHPC Behavior
Five parameters relating to plasticity response in UHPC, namely σ b0 /σ c0 , k c , ψ, ξ, µ, are required to define CDP model. Two parameters of σ b0 /σ c0 and k c modify the yield surface. σ b0 /σ c0 is the ratio of biaxial compressive strength to uniaxial compressive strength which influences the yield surface in a plane stress state and the parameter k c is the ratio between distances measured from the hydrostatic axis to tensile and compressive meridians and is used to define the shape of the failure surface in deviatoric plane. The other two parameters, ψ and ξ modify the non-associated potential flow. ψ is the dilation angle that describes the angle of inclination of the failure surface toward the hydrostatic axis measured in the meridional plane. ξ is an eccentricity parameter that controls the deviation of the hyperbolic plastic potential from its asymptote. µ is the viscosity parameter which is used for the visco-plastic regularization of the concrete constitutive equations [33,35,39,43]. These five parameters were selected to be 1.05, 2/3, 30, 0.1, 1E-4, respectively, according to the literature and sensitivity analysis conducted in this study [33,35,38,[43][44][45]. In addition, Poisson's ratio of UHPC (or UHPFRC) were considered to be 0.2 according to AFGC [8]. Density of UHPC (or UHPFRC) was measured to be 2565 kg/m 3 [17].
CDP model accounts for stiffness degradation in compression and tension. This stiffness degradation is defined in terms of scalar variables that are functions of the plastic strains. These variables capture degradation in material stiffness with increased loading. In other words, the degradation variables get more pronounced with increasing plastic strain which are zero and unity (=1) for an undamaged state and a complete damage state, respectively.
The tensile damage parameter (d t ) gets activated after attaining peak tensile strength. Therefore, damage contours indicate tensile cracking and the level of damage increases with increasing strain (crack widening) at higher load levels. To account for degradation in stiffness due to cracking, tensile damage parameter (shown in Equation (1)) was included in the model. Moreover, the stiffness reduction in compression was incorporated in the plasticity model as presented in Equation (2) [34,35,43].
where in σ t , f t , E, σ c , d t , and d c are tensile stress, tensile strength, elastic modulus, compressive stress, tensile, and compressive damage parameters, respectively.

Compressive Behavior of UHPC
There are limited test data and associated relations for tracing uniaxial compressive stress-strain behavior of UHPC and UHPFRC. A material model that can relate stressstrain behavior of UHPC (or UHPFRC) to compressive strength and modulus of elasticity is developed for undertaking detailed analysis of UHPC and UHPFRC structures.
The uniaxial stress-strain response of UHPFRC can be approximated with a quad-part model that include softening branch as well. A schematic representation of the proposed model of stress-strain response of UHPFRC under compression is shown in Figure 2a. UHPFRC under compression, unlike conventional concrete, exhibits nearly linear behavior up to almost 70% of their compressive strength [46] (point 2 in Figure 2a). Compressive strength of UHPFRC which varies depending on fiber content and type, curing regime, mix design, etc., is to be determined by uniaxial compression test. The stress-strain response in this linear part can be defined by elastic modulus using Equation (3). The linear part of stress strain response is followed by a nonlinear phase until peak strength is reached and can be represented using Equation (4), which relates strain to stress by elastic modulus and a reduction factor (α). This reduction factor (Equation (5)), defines reduction of stress from linear elastic stress [46].
in which α the reduction factor is given as: where ε, E, f c , and α are compressive strain, elastic modulus, compressive strength, and reduction factor, respectively. Using data from literature [13,34,[46][47][48][49][50][51][52][53][54], elastic modulus of UHPFRC is plotted as a function of ( in Figure 2b and an empirical relation available in literature for calculating elastic modulus of high strength concrete [55] is modified for UHPFRC as: curing regime, mix design, etc., is to be determined by uniaxial compression test. The stress-strain response in this linear part can be defined by elastic modulus using Equation (3). The linear part of stress strain response is followed by a nonlinear phase until peak strength is reached and can be represented using Equation (4), which relates strain to stress by elastic modulus and a reduction factor (α). This reduction factor (Equation (5)), defines reduction of stress from linear elastic stress [46]. The ascending branch of compressive stress strain response (1-2, and 2-3) is calculated using Equations (3)- (6). The descending branch (3-4, and 4-5), is obtained based on empirically derived parameters generated in experiments by Empelmann [3]. The relations to calculate these five key points based on fiber content (by volume) and fiber size (aspect ratio) are given in Figure 2a. Where ε 0 is the strain corresponding to compressive strength and v f , l f , and d f are fiber volume fraction, fiber length, and diameter, respectively. Behavior of UHPC with no fibers is also linear up to 70% of compressive strength and it fails in brittle manner under compression (explosive) [56] as shown in Figure 2a. Predicted stress-strain response of UHPFRC based on parameters generated from two different experiments is plotted along with experimental stress-strain response in Figure 3. It shows a good agreement between the predicted stress-strain response using proposed model and test data especially in ascending branch of response.
Modelling 2021, 2, FOR PEER REVIEW 5 = (1 − ) 0.70 ′ < < ′ (4) in which α the reduction factor is given as: where , , ′ , and are compressive strain, elastic modulus, compressive strength, and reduction factor, respectively. Using data from literature [13,34,[46][47][48][49][50][51][52][53][54], elastic modulus of UHPFRC is plotted as a function of ( The ascending branch of compressive stress strain response (1-2, and 2-3) is calculated using Equations (3)- (6). The descending branch (3-4, and 4-5), is obtained based on empirically derived parameters generated in experiments by Empelmann [3]. The relations to calculate these five key points based on fiber content (by volume) and fiber size (aspect ratio) are given in Figure 2a. Where 0 is the strain corresponding to compressive strength and , , and are fiber volume fraction, fiber length, and diameter, respectively. Behavior of UHPC with no fibers is also linear up to 70% of compressive strength and it fails in brittle manner under compression (explosive) [56] as shown in Figure 2a. Predicted stress-strain response of UHPFRC based on parameters generated from two different experiments is plotted along with experimental stress-strain response in Figure 3. It shows a good agreement between the predicted stress-strain response using proposed model and test data especially in ascending branch of response.

Tensile Behavior of UHPC
To capture the beneficial effects of UHPFRC, tensile response of UHPFRC is to be properly replicated before and after cracking. For this, stress-strain response in uniaxial tension is needed to capture the hardening and softening response of concrete.
Behavior of UHPC without fibers under tension after cracking is brittle and does not exhibit strain hardening and a significant descending branch [56] as shown in Figure  4a. However, the fibers present in UHPFRC induce significant bridging stress between open crack faces. This bridging stress between opened cracks faces leads to a relatively higher fracture toughness and ductility in UHPFRC. Therefore, it is essential that this fi-

Tensile Behavior of UHPC
To capture the beneficial effects of UHPFRC, tensile response of UHPFRC is to be properly replicated before and after cracking. For this, stress-strain response in uniaxial tension is needed to capture the hardening and softening response of concrete.
Behavior of UHPC without fibers under tension after cracking is brittle and does not exhibit strain hardening and a significant descending branch [56] as shown in Figure 4a. However, the fibers present in UHPFRC induce significant bridging stress between open crack faces. This bridging stress between opened cracks faces leads to a relatively higher fracture toughness and ductility in UHPFRC. Therefore, it is essential that this fiber bridging mechanism is effectively incorporated in modeling tensile fracture of UHPFRC through stress-strain response in tension. Typical stress strain behavior of UHPFRC is idealized into three stages as shown in Figure 4a [57]. The initial part is linear elastic up to cracking stress of UHPFRC, which is followed by strain hardening part accompanied by initiation of multiple cracking facilitated by fiber bridging. This is further followed by softening branch that presents crack opening with fiber bridging.
Modelling 2021, 2, FOR PEER REVIEW 6 FRC through stress-strain response in tension. Typical stress strain behavior of UHPFRC is idealized into three stages as shown in Figure 4a [57]. The initial part is linear elastic up to cracking stress of UHPFRC, which is followed by strain hardening part accompanied by initiation of multiple cracking facilitated by fiber bridging. This is further followed by softening branch that presents crack opening with fiber bridging. Tensile behavior of UHPFRC cast for this study [17] as well as UHPFRC studied by Sing et al. [34] and Wille et al. [57] under direct tension is presented in Figure 4b. It can be seen that there is significant variation in tensile stress-strain response of UHPFRC with different mix, fiber type, volume fraction, and distribution. Since the shape of tensile stress-strain response of UHPFRC is highly influenced by various factors such as characteristic strength of the concrete matrix, fiber type, orientation and distribution of fibers, fiber aspect-ratio and fiber content (by volume), the key points on stress strain response plot, namely cracking, peak, and ultimate points are to be evaluated from direct tension tests. In addition, uniaxial tensile strength can be calculated through flexural tests on concrete prisms. However, flexural capacity alone cannot be utilized to describe the complete fracture mechanism. Therefore, inverse analysis is to be conducted to develop a tensile fracture model for UHPFRC from flexural tests [1,2,49,51,58].

Steel Reinforcement
An available plasticity model in ABAQUS based on Mises yield surface with associated plastic flow and isotropic hardening was utilized for constitutive modeling of steel reinforcing bars. The response of steel bars under tension and compression consisting of three phases of linear elastic, yield plateau, and strain hardening is incorporated in the model. The strain hardening part of stress strain curve can be calculated using Equation (7) [59].
where , , ℎ , and are yield strength, ultimate strength, strain at the end of the plateau part, and ultimate strain, respectively. Tensile behavior of UHPFRC cast for this study [17] as well as UHPFRC studied by Sing et al. [34] and Wille et al. [57] under direct tension is presented in Figure 4b. It can be seen that there is significant variation in tensile stress-strain response of UHPFRC with different mix, fiber type, volume fraction, and distribution. Since the shape of tensile stressstrain response of UHPFRC is highly influenced by various factors such as characteristic strength of the concrete matrix, fiber type, orientation and distribution of fibers, fiber aspect-ratio and fiber content (by volume), the key points on stress strain response plot, namely cracking, peak, and ultimate points are to be evaluated from direct tension tests. In addition, uniaxial tensile strength can be calculated through flexural tests on concrete prisms. However, flexural capacity alone cannot be utilized to describe the complete fracture mechanism. Therefore, inverse analysis is to be conducted to develop a tensile fracture model for UHPFRC from flexural tests [1,2,49,51,58].

Steel Reinforcement
An available plasticity model in ABAQUS based on Mises yield surface with associated plastic flow and isotropic hardening was utilized for constitutive modeling of steel reinforcing bars. The response of steel bars under tension and compression consisting of three phases of linear elastic, yield plateau, and strain hardening is incorporated in the model. The strain hardening part of stress strain curve can be calculated using Equation (7) [59].
where f y , f u , ε sh , and ε u are yield strength, ultimate strength, strain at the end of the plateau part, and ultimate strain, respectively.

Bond Slip Behavior of UHPC
Number of researchers have investigated bond behavior of steel reinforcement with UHPC. These authors reported that the average bond strength of steel bars embedded in UHPC is about ten times higher than that of steel bars placed in conventional concrete [50].
The modified CMR [60] model proposed by Yoo et al. [50] for UHPC and UHPFRC with varying fiber content was utilized to define bond behavior of steel bars embedded in UHPC. The CMR model is defined in Equations (8) and (9), where τ and τ max are bond stress and bond strength, respectively. s r and β are the coefficients based on curve fitting of experimental results on UHPC and UHPFRC which were adopted as 0.07, and 0.8, respectively.
The bond force between the concrete and reinforcing steel bar for the bond element is obtained by multiplying contact area between reinforcing bar and concrete by average bond stress between concrete and steel reinforcement.

Analysis Details
Since UHPFRC beams, owing to high ductility, undergo large deflections (as compared to beams made of conventional concrete or UHPC), the effect of geometric non-linearity should be given due consideration in the analysis. This was considered in the model through updated Lagrangian method [39]. Newton-Raphson solution technique is adopted and a tolerance limit of 0.5% of average force is specified to meet convergence criterion at each load increment [61]. The analysis is carried out in small displacement increments which are automatically chosen by ABAQUS. The mesh size was arrived by carrying out a parametric study with different mesh sizes and a mesh size of 25 mm was able to predict good post yield response in the beams selected for validation.

Modeling Interfacial Bond between Rebar and Concrete
The interfacial bond between reinforcement and concrete can be accounted for through explicit modeling of both reinforcement ribs and the concrete lugs [62]. Alternatively, local bond-slip can be modeled as bond-link elements which provides a reasonable compromise between accuracy and computational efficiency [40].
In bond-link element approach, the concrete and the reinforcing steel are represented by two different sets of elements, and node pairs at the interface are connected using interfacial spring elements. Three spring elements are modeled at each node pair, wherein one spring represents shear bond behavior according to a bond-slip relationship. This bond-slip is related to longitudinal axis direction. The other two springs represent the bond behavior in the normal directions which are assumed to be rigid by assigning large spring stiffness to the normal springs [40].

Model Validation
The developed finite element-based model was validated by analyzing a set of UHPC and UHPFRC beams tested in the laboratory [5,11,13,17,34]. To evaluate the effectiveness of the developed model in tracing structural response of UHPC and UHPFRC beams, model predictions including load carrying capacity, load-strain behavior, load-deflection behavior, and cracking were compared with results and observations from experiments [63].

Selected Beams for Validation
Four UHPFRC beams, designated as U-B3 to U-B6 were fabricated and tested for tracing flexural and shear response of beams. Two beams (U-B3 and U-B5) were tested under flexure by applying two point loads on top face of the beams as shown in Figure 5. This test set-up ensured that the critical span between loading points was subjected to pure flexure and no shear. Two beams (U-B4 and U-B6) were subjected to dominant shear loading through applying a single point load at the distance of 610 mm from support (see Figure 5) [17,64].
Modelling 2021, 2, FOR PEER REVIEW 8 shear loading through applying a single point load at the distance of 610 mm from support (see Figure 5) [17,64].

Figure 5.
Test set-up, and cross sectional details of tested UHPFRC beams [17].
In addition, results from tests on UHPC and UHPFRC beams conducted by Yoo and Yoon [11], Yang et al. [5], and Singh et al. [34] were utilized to validate the model. Details of these beam are presented in Table 1. Additional details on experiments, including arrangement of reinforcing bars, loading set-up, and material properties can be found in relevant references [7,11,13,17,34,64]. The compressive strength of UHPC and UHPFRC cylinders reported by authors were used to define compressive stress strain response as proposed in Section 2.2.2. Properties of concrete in tension were based on results of the tests conducted under direct tension. Stress strain response of reinforcement bars were calculated using Equation (7) based on reported yield and ultimate stress of the reinforcement.  (1) and (2) are peak loads obtained from FEA and tests, respectively.

Load-Deflection Response
To establish the validity of the model, a number of UHPC and UHPFRC beams, as listed in Table 1, with different material characteristics, and tested under different load- Figure 5. Test set-up, and cross sectional details of tested UHPFRC beams (data from [17]).
In addition, results from tests on UHPC and UHPFRC beams conducted by Yoo and Yoon [11], Yang et al. [5], and Singh et al. [34] were utilized to validate the model. Details of these beam are presented in Table 1. Additional details on experiments, including arrangement of reinforcing bars, loading set-up, and material properties can be found in relevant references [7,11,13,17,34,64]. The compressive strength of UHPC and UHPFRC cylinders reported by authors were used to define compressive stress strain response as proposed in Section 2.2.2. Properties of concrete in tension were based on results of the tests conducted under direct tension. Stress strain response of reinforcement bars were calculated using Equation (7) based on reported yield and ultimate stress of the reinforcement.

Load-Deflection Response
To establish the validity of the model, a number of UHPC and UHPFRC beams, as listed in Table 1, with different material characteristics, and tested under different loading conditions, were analyzed and load-deflection response of the beams are compared with measured data from experiments. Figure 6a,b shows the load-deflection response of beams U-B4 and U-B5 [17] with different reinforcement ratios tested under different loading configurations (dominant shear and flexure) while Figure 6c,d shows load deflection of beams ρ 0.94%-NF and B25 [11,34] with different fiber volume fractions (0, and 2.25%). It can be seen from the overall trends that the prediction from the numerical model compares well with the measured data in different stages i.e., elastic stage till tensile cracking initiation, post-cracking stage with propagation of cracks, yielding in steel reinforcing bars, and plastic deformation stage including softening until failure in both cases of UHPC and UHPFRC beams. It should be noted that similar comparisons were made for number of other beams tested by various researchers and very good comparisons were obtained (see Table 1).
ing conditions, were analyzed and load-deflection response of the beams are compared with measured data from experiments. Figure 6a,b shows the load-deflection response of beams U-B4 and U-B5 [17] with different reinforcement ratios tested under different loading configurations (dominant shear and flexure) while Figure 6c,d shows load deflection of beams 0.94%-NF and B25 [11,34] with different fiber volume fractions (0, and 2.25%). It can be seen from the overall trends that the prediction from the numerical model compares well with the measured data in different stages i.e., elastic stage till tensile cracking initiation, post-cracking stage with propagation of cracks, yielding in steel reinforcing bars, and plastic deformation stage including softening until failure in both cases of UHPC and UHPFRC beams. It should be noted that similar comparisons were made for number of other beams tested by various researchers and very good comparisons were obtained (see Table 1).
Moreover, the developed model predicts response well for both UHPC and UHP-FRC beams under different loading conditions i.e., flexural and dominant shear loading. Only, in post cracking stage, model predictions are slightly stiffer than the measured experimental data. This difference can be attributed to possible cracks developing in concrete due to dry shrinkage and environmental effects which resulted in softer response in tests. UHPFRC beams with higher reinforcement ratio as compared to similar beams with lower reinforcement ratio (beams U-B5 and U-B6 as compared to beams U-B3 and U-B4), exhibited higher strain hardening in experiments as compared to finite element model predictions. This can be attributed to slight variations in the stress-strain curve adopted for reinforcing steel, which gets more pronounced in load-deflection response  Moreover, the developed model predicts response well for both UHPC and UHPFRC beams under different loading conditions i.e., flexural and dominant shear loading. Only, in post cracking stage, model predictions are slightly stiffer than the measured experimental data. This difference can be attributed to possible cracks developing in concrete due to dry shrinkage and environmental effects which resulted in softer response in tests.
UHPFRC beams with higher reinforcement ratio as compared to similar beams with lower reinforcement ratio (beams U-B5 and U-B6 as compared to beams U-B3 and U-B4), exhibited higher strain hardening in experiments as compared to finite element model predictions. This can be attributed to slight variations in the stress-strain curve adopted for reinforcing steel, which gets more pronounced in load-deflection response at higher reinforcement ratios.
In addition to load-deflection response, the failure load of the beams, as predicted from FEA, is compared with the measured peak loads from experiments in Table 1. A ratio of unity (=1) indicates perfect agreement, while less than unity and higher than unity are conservative and unconservative predictions, respectively. Ratio of total load carrying capacity (P) from model to that of measured values in tests ranges from 0.90 to 1.03 for analyzed UHPC and UHPFRC beams, indicating good capability of the model to capture failure load.
The model is also capable of evaluating the bond developed in rebars at different stages of loading. The bond developed between steel reinforcement and concrete in beams U-B5 and U-B6 (subjected to flexural and dominant shear loading) is plotted in Figure 7. It can be seen that the bond stress increased from about 0.4 MPa, just prior to cracking initiated in the beams, to a stress of about 5.5 MPa at failure. The maximum bond stress in beams U-B5 and U-B6 developed at critical section of the beams i.e., under load points.
odelling 2021, 2, FOR PEER REVIEW 1 In addition to load-deflection response, the failure load of the beams, as predicted from FEA, is compared with the measured peak loads from experiments in Table 1. A ra tio of unity (=1) indicates perfect agreement, while less than unity and higher than unity are conservative and unconservative predictions, respectively. Ratio of total load carry ing capacity (P) from model to that of measured values in tests ranges from 0.90 to 1.03 for analyzed UHPC and UHPFRC beams, indicating good capability of the model to capture failure load.
The model is also capable of evaluating the bond developed in rebars at differen stages of loading. The bond developed between steel reinforcement and concrete in beams U-B5 and U-B6 (subjected to flexural and dominant shear loading) is plotted in Figure 7. It can be seen that the bond stress increased from about 0.4 MPa, just prior to cracking initiated in the beams, to a stress of about 5.5 MPa at failure. The maximum bond stress in beams U-B5 and U-B6 developed at critical section of the beams i.e., unde load points.

Load-Strain Response
To illustrate the capability of model in capturing the local behavior of UHPFRC beams, predicted load-longitudinal strains on reinforcing bar and concrete layers at criti cal section of beams U-B3 and U-B4 [17] are compared with the measured strains from tests in Figure 8. A negative strain in the figure indicates that the material is in compres sive state while positive strain indicates the presence of tensile state. Strain predictions are plotted until strain gauges stopped recording strains reliably due to cracks develop ing at the location of the strain gauges. The trends in the plots show that predicted strains from the model are in good agreement with the measured strains from tests However, predicted strain response in post cracking state is stiffer as compared to the measured experimental values. This difference between predictions from model and measured data from experiments can be due to variations resulting from difference in the material models, as well as any experimental discrepancies arising from bond leve between concrete and strain gauges.

Load-Strain Response
To illustrate the capability of model in capturing the local behavior of UHPFRC beams, predicted load-longitudinal strains on reinforcing bar and concrete layers at critical section of beams U-B3 and U-B4 [17] are compared with the measured strains from tests in Figure 8. A negative strain in the figure indicates that the material is in compressive state while positive strain indicates the presence of tensile state. Strain predictions are plotted until strain gauges stopped recording strains reliably due to cracks developing at the location of the strain gauges. The trends in the plots show that predicted strains from the model are in good agreement with the measured strains from tests. However, predicted strain response in post cracking state is stiffer as compared to the measured experimental values. This difference between predictions from model and measured data from experiments can be due to variations resulting from difference in the material models, as well as any experimental discrepancies arising from bond level between concrete and strain gauges.

Crack Propagation and Failure Mode
The developed numerical model is also capable of capturing crack progression through tracing scalar tensile damage parameter. This damage parameter in tension gets activated after concrete attains its peak tensile stress. Therefore, damage contours indicate tensile cracking and the level of damage increases with increasing strains (crack widening) at higher load levels. In other words, tensile damage parameter of 0 (zero) represents no tension damage, while a value of 1 (unity) represents complete damage state. Direction of cracking can be represented using direction of principal strains being perpendicular to crack direction.
Tensile damage contours and crack direction obtained through FEA along with experimental results are shown in Figure 9 for beams U-B5 and U-B6, subjected to differen loading configurations. In beam U-B5, subjected to flexural loading, flexural cracks initiated at load level of 28 kN at extreme tension fibers of the beam. These flexural cracks developed at the region between loading points being under bending. At a load level o 86 kN, this tensile damage propagated toward compression zone (till mid depth of the beam) as shown in Figure 9a. Direction of principal strain generated from the mode confirms propagation of flexural cracks. This type of cracking behavior was observed in the experiment, also can be seen in Figure 9a. With increasing load, at P = 108 kN which is close to failure load, greater depth of the beam was subjected to tensile damage. This clearly infers propagation of cracking toward compression zone with increasing load as observed in the tests (Figure 9a). The tension damage contour generated through FEA illustrates that maximum tensile damage, close to failure, is concentrated at the critica section of the beam which matches with the progression of macro crack as observed during the experiment. It should be noted that a greater value of tensile damage parameter represents a greater level of damage and larger crack width.
In the case of UHPFRC beam U-B6, subjected to shear dominant loading, tensile damage got initiated at extreme tension fiber at a load level of 40 kN and this resulted from tensile stress exceeding the tensile strength of concrete. This initial tensile damage (cracking) was confined to the zone between the load point and mid-span and was mainly in the form of flexural cracks in lower depth of the beam. When the load on the beam increased to 137 kN, these flexural cracks propagated toward compression zone (upper depth). Further, additional flexural-shear cracks developed in the shear span as shown in Figure 9b. The tensile damage contour and principal direction shown in Figure   Figure 8. Comparison of load-strain response in bars and concrete at critical section of beams, as obtained from FEA and tests.

Crack Propagation and Failure Mode
The developed numerical model is also capable of capturing crack progression through tracing scalar tensile damage parameter. This damage parameter in tension gets activated after concrete attains its peak tensile stress. Therefore, damage contours indicate tensile cracking and the level of damage increases with increasing strains (crack widening) at higher load levels. In other words, tensile damage parameter of 0 (zero) represents no tension damage, while a value of 1 (unity) represents complete damage state. Direction of cracking can be represented using direction of principal strains being perpendicular to crack direction.
Tensile damage contours and crack direction obtained through FEA along with experimental results are shown in Figure 9 for beams U-B5 and U-B6, subjected to different loading configurations. In beam U-B5, subjected to flexural loading, flexural cracks initiated at load level of 28 kN at extreme tension fibers of the beam. These flexural cracks developed at the region between loading points being under bending. At a load level of 86 kN, this tensile damage propagated toward compression zone (till mid depth of the beam) as shown in Figure 9a. Direction of principal strain generated from the model confirms propagation of flexural cracks. This type of cracking behavior was observed in the experiment, also can be seen in Figure 9a. With increasing load, at P = 108 kN which is close to failure load, greater depth of the beam was subjected to tensile damage. This clearly infers propagation of cracking toward compression zone with increasing load as observed in the tests (Figure 9a). The tension damage contour generated through FEA illustrates that maximum tensile damage, close to failure, is concentrated at the critical section of the beam which matches with the progression of macro crack as observed during the experiment. It should be noted that a greater value of tensile damage parameter represents a greater level of damage and larger crack width. ed further toward compression zone and more cracks got initiated and propagated further into the shear span. The predicted tensile damage and principal direction in the shear span of the beam U-B6 (between left support and loading point) close to failure and cracking pattern observed in experiment are compared in Figure 9b. The predicted direction of principal strains in shear span matches well with the direction of the major diagonal tension crack (which is perpendicular to principal direction) observed in the experiment. Figure 9. Comparison of tensile damage contours and principal strain direction as predicted by FEA with test data.

Concrete and Stirrups Contribution to Shear Capacity
UHPFRC possesses high tensile strength and ductile characteristics (high ultimate tensile strain), and this can be utilized to realize high shear capacity in UHPFRC beams. The developed model was applied to explore the feasibility of removing shear reinforcement in UHPFRC beams; specifically, the contribution of concrete and stirrups to shear resistance was quantified. Behavior of tested UHPFRC beams U-B4 and U-B6, provided with shear reinforcement, were analyzed under dominant shear loading using the above developed numerical model. Moreover, for comparative study, NSC beams with the same cross sectional details as U-B4 and U-B6, as shown in Figure 5, were modeled under dominant shear loading. The material model for NSC recommended in ABAQUS documentation [39] were incorporated into the model for analyzing NSC beams. In the case of UHPFRC beam U-B6, subjected to shear dominant loading, tensile damage got initiated at extreme tension fiber at a load level of 40 kN and this resulted from tensile stress exceeding the tensile strength of concrete. This initial tensile damage (cracking) was confined to the zone between the load point and mid-span and was mainly in the form of flexural cracks in lower depth of the beam. When the load on the beam increased to 137 kN, these flexural cracks propagated toward compression zone (upper depth). Further, additional flexural-shear cracks developed in the shear span as shown in Figure 9b. The tensile damage contour and principal direction shown in Figure 9b match well with the crack pattern observed during tests. As the load increased further to 155 kN (just prior to failure), shear stresses in the shear span increased significantly and shear cracks became much more predominant. In other words, maximum principal stresses in shear span exceeded the tensile capacity of UHPFRC. Model predictions also show that as the load level approached failure load, tensile damage propagated further toward compression zone and more cracks got initiated and propagated further into the shear span. The predicted tensile damage and principal direction in the shear span of the beam U-B6 (between left support and loading point) close to failure and cracking pattern observed in experiment are compared in Figure 9b. The predicted direction of principal strains in shear span matches well with the direction of the major diagonal tension crack (which is perpendicular to principal direction) observed in the experiment.

Concrete and Stirrups Contribution to Shear Capacity
UHPFRC possesses high tensile strength and ductile characteristics (high ultimate tensile strain), and this can be utilized to realize high shear capacity in UHPFRC beams. The developed model was applied to explore the feasibility of removing shear reinforcement in UHPFRC beams; specifically, the contribution of concrete and stirrups to shear resistance was quantified. Behavior of tested UHPFRC beams U-B4 and U-B6, provided with shear reinforcement, were analyzed under dominant shear loading using the above developed numerical model. Moreover, for comparative study, NSC beams with the same cross sectional details as U-B4 and U-B6, as shown in Figure 5, were modeled under dominant shear loading. The material model for NSC recommended in ABAQUS documentation [39] were incorporated into the model for analyzing NSC beams.
Predicted load-deflection response of NSC beams, with different longitudinal tensile reinforcement ratios (i.e., ρ t = 0.90% and ρ t = 1.20%) with and without stirrups, is shown in Figure 10a,b. As can be seen in the figure, the load carrying capacity in the NSC beams, with tensile reinforcement ratio of 0.90% and 1.20%, reduces by 14% and 15% when the stirrups are removed. In addition, NSC beams without shear reinforcement exhibited significant reduction of stiffness after reaching peak load, as compared to the beams provided with stirrups.
Modelling 2021, 2, FOR PEER REVIEW 13 Predicted load-deflection response of NSC beams, with different longitudinal tensile reinforcement ratios (i.e., = 0.90% and = 1.20%) with and without stirrups, is shown in Figure 10a,b. As can be seen in the figure, the load carrying capacity in the NSC beams, with tensile reinforcement ratio of 0.90% and 1.20%, reduces by 14% and 15% when the stirrups are removed. In addition, NSC beams without shear reinforcement exhibited significant reduction of stiffness after reaching peak load, as compared to the beams provided with stirrups. Moreover, stress distribution developed along the tensile reinforcement bars in NSC beams, with and without stirrups, at peak load level is shown in Figure 11a,b. The comparative trends indicate that reinforcing bars in NSC beam provided with stirrups, yielded as opposed to NSC beam with no stirrups; in other words, the beam with stirrups reached its ultimate moment capacity. However, NSC beam without stirrups failed in shear before reaching ultimate moment capacity (see Figure 11a,b). Moreover, stress distribution developed along the tensile reinforcement bars in NSC beams, with and without stirrups, at peak load level is shown in Figure 11a,b. The comparative trends indicate that reinforcing bars in NSC beam provided with stirrups, yielded as opposed to NSC beam with no stirrups; in other words, the beam with stirrups reached its ultimate moment capacity. However, NSC beam without stirrups failed in shear before reaching ultimate moment capacity (see Figure 11a,b). The behavior of UHPFRC beams having same configurations as that of U-B4 and U-B6 was analyzed under two cases; one with stirrups and the other one without stirrups. The load-deflection response of UHPFRC beams with and without shear reinforcement are compared in Figure 10c,d. It can be seen that removing stirrups did not affect the overall structural response of the beams in terms of load-carrying capacity and ductility. This was also confirmed with experimental data generated from tests on beams U-B3 to U-B6. Beams U-B4 and U-B6, with no shear reinforcement reached their ultimate moment capacity under dominant shear loading.
The stress developed in longitudinal bars in UHPFRC beam (similar to U-B4) with and without stirrups is presented in Figure 11c,d. In both cases of UHPFRC beam with and without stirrups, bars yielded. In other words, unlike NSC beams, UHPFRC beams without stirrups do not experience abrupt failure before yielding in the reinforcing bars and this is owing to high tensile strength and ultimate strain of UHPFRC that develops due to bridging effect facilitated by presence of steel fibers.
The contribution of different components to shear capacity of NSC and UHPFRC beams was quantified to determine the extent of contribution from UHPFRC and stirrups. For this purpose, the stirrups crossing crack surface, which were in tension as compared to other stirrups, were identified. The contribution of stirrups to shear strength (Vs) was calculated using the tensile stress developed in the stirrups of two NSC and UHPFRC beams. In the case of a NSC beam contribution of stirrups to shear resistance, which was small till cracking, it increased after first cracking and reached to about 90% of total shear capacity at peak state (see Figure 12). However, in the case of a The behavior of UHPFRC beams having same configurations as that of U-B4 and U-B6 was analyzed under two cases; one with stirrups and the other one without stirrups. The load-deflection response of UHPFRC beams with and without shear reinforcement are compared in Figure 10c,d. It can be seen that removing stirrups did not affect the overall structural response of the beams in terms of load-carrying capacity and ductility. This was also confirmed with experimental data generated from tests on beams U-B3 to U-B6. Beams U-B4 and U-B6, with no shear reinforcement reached their ultimate moment capacity under dominant shear loading.
The stress developed in longitudinal bars in UHPFRC beam (similar to U-B4) with and without stirrups is presented in Figure 11c,d. In both cases of UHPFRC beam with and without stirrups, bars yielded. In other words, unlike NSC beams, UHPFRC beams without stirrups do not experience abrupt failure before yielding in the reinforcing bars and this is owing to high tensile strength and ultimate strain of UHPFRC that develops due to bridging effect facilitated by presence of steel fibers.
The contribution of different components to shear capacity of NSC and UHPFRC beams was quantified to determine the extent of contribution from UHPFRC and stirrups. For this purpose, the stirrups crossing crack surface, which were in tension as compared to other stirrups, were identified. The contribution of stirrups to shear strength (Vs) was calculated using the tensile stress developed in the stirrups of two NSC and UHPFRC beams. In the case of a NSC beam contribution of stirrups to shear resistance, which was small till cracking, it increased after first cracking and reached to about 90% of total shear capacity at peak state (see Figure 12). However, in the case of a UHPFRC beam, shear resistance of the beam resulted mostly from contribution of concrete (Vc) and stirrups did not contribute to shear capacity of the beam as shown in Figure 12.
Modelling 2021, 2, FOR PEER REVIEW 15 UHPFRC beam, shear resistance of the beam resulted mostly from contribution of concrete (Vc) and stirrups did not contribute to shear capacity of the beam as shown in Fig Contribution of concrete to shear strength comprises contribution of uncracked concrete in compression (from compression block of concrete), and resistance arising from fiber bridging and aggregate interlock across cracked concrete. It should be noted that in the case of slender steel fiber-reinforced concrete beams without stirrups, and under reinforced beams, dowel action contribution can be neglected [65].
Therefore, to quantify the contribution of uncracked concrete (Vcc) to shear resistance, model prediction is used to identify the region above the neutral axis which is in compression (see Figure 13). This region (compression block) is subjected to normal compressive and shear stresses. Therefore, the contribution of concrete to shear strength arising from compression block can be evaluated utilizing predicted shear stress in region above neutral axis. Once the contribution of uncracked concrete to shear capacity is determined, the remaining shear strength of concrete is attributed to interfacial shear resistance (Vi) which arises from aggregate interlock and fiber bridging as shown in Figure  13 [65].
In two UHPFRC beams without stirrups, U-B4 and U-B6, subjected to dominant shear loading, shear strength contribution from compression block (Vcc) and fiber bridging and aggregate interlock (Vi) is evaluated and presented in Table 2. The results show that 67% and 65% of shear capacity of these beams (U-B4 and U-B6) came from compression block (Vcc) at initial cracking stage. With increasing load level (at peak and failure load levels), cracks propagated more toward the compression zone. Therefore, the con- Contribution of concrete to shear strength comprises contribution of uncracked concrete in compression (from compression block of concrete), and resistance arising from fiber bridging and aggregate interlock across cracked concrete. It should be noted that in the case of slender steel fiber-reinforced concrete beams without stirrups, and under reinforced beams, dowel action contribution can be neglected [65].
Therefore, to quantify the contribution of uncracked concrete (Vcc) to shear resistance, model prediction is used to identify the region above the neutral axis which is in compression (see Figure 13). This region (compression block) is subjected to normal compressive and shear stresses. Therefore, the contribution of concrete to shear strength arising from compression block can be evaluated utilizing predicted shear stress in region above neutral axis. Once the contribution of uncracked concrete to shear capacity is determined, the remaining shear strength of concrete is attributed to interfacial shear resistance (Vi) which arises from aggregate interlock and fiber bridging as shown in Figure 13 [65].
In addition, contribution of compression block and aggregate interlock to shear strength of NSC beams, having the same cross sectional details as beams U-B4 and U-B6 is quantified (see Table 2). The results show that NSC beams did not exhibit significant reduction in compression block contribution to shear strength (Vcc/V) with increasing load till peak state. In other words, contribution of interfacial shear resistance to shear capacity (Vi/V) did not increase, as opposed to UHPFRC beams wherein interfacial shear resistance (Vi/V) increased due to activation of fiber bridging. Figure 13. Schematic of assumed strain distribution and internal stresses in shear span of a beam [65].
These analyses illustrate the usefulness of the developed model in determining the contribution of concrete and stirrups to shear capacity of beams. In UHPFRC beams, due to high tensile and compression strength of UHPFRC, concrete contribution to shear capacity is quite significant and until final stages of loading, stirrups do not play major role in resisting shear. Therefore, stirrups can be eliminated in many cases. This is in contrast to the shear response of NSC beams wherein contribution of concrete to shear capacity decreases after the onset of cracks in beams, and thus the presence of stirrups are critical in resisting shear beyond cracking load levels. Figure 13. Schematic of assumed strain distribution and internal stresses in shear span of a beam (adapted from [65]).
In two UHPFRC beams without stirrups, U-B4 and U-B6, subjected to dominant shear loading, shear strength contribution from compression block (Vcc) and fiber bridging and aggregate interlock (Vi) is evaluated and presented in Table 2. The results show that 67% and 65% of shear capacity of these beams (U-B4 and U-B6) came from compression block (Vcc) at initial cracking stage. With increasing load level (at peak and failure load levels), cracks propagated more toward the compression zone. Therefore, the contribution of compression block to shear capacity decreased as smaller depth of concrete was in compression. As can be seen in Table 2, with decreasing contribution of compression block (Vcc/V), interfacial shear resistance (Vi/V) which is due to fiber bridging and aggregate interlock at cracks surfaces increased. The reduction in contribution of compression block to shear resistance was higher in case of beam U-B6 as compared to beam U-B4. This is attributed to higher applied load in beam U-B6 as compared to beam U-B4 due to higher reinforcement ratio. Therefore, cracks propagated more toward compression zone in beam U-B6, and smaller depth of the beam was uncracked (in compression). In addition, contribution of compression block and aggregate interlock to shear strength of NSC beams, having the same cross sectional details as beams U-B4 and U-B6 is quantified (see Table 2). The results show that NSC beams did not exhibit significant reduction in compression block contribution to shear strength (Vcc/V) with increasing load till peak state. In other words, contribution of interfacial shear resistance to shear capacity (Vi/V) did not increase, as opposed to UHPFRC beams wherein interfacial shear resistance (Vi/V) increased due to activation of fiber bridging.
These analyses illustrate the usefulness of the developed model in determining the contribution of concrete and stirrups to shear capacity of beams. In UHPFRC beams, due to high tensile and compression strength of UHPFRC, concrete contribution to shear capacity is quite significant and until final stages of loading, stirrups do not play major role in resisting shear. Therefore, stirrups can be eliminated in many cases. This is in contrast to the shear response of NSC beams wherein contribution of concrete to shear capacity decreases after the onset of cracks in beams, and thus the presence of stirrups are critical in resisting shear beyond cracking load levels.

Conclusions
Based on the results presented in the paper the following conclusions can be drawn.
• UHPC or UHPFRC exhibits significantly different mechanical properties as compared to conventional concrete. The proposed numerical model utilizing concrete damage plasticity model with adjusted parameters is capable of tracing the structural response of UHPC beams in the entire range of loading; from precracking stage till failure. • The numerical model presented here can accommodate various configurations in beams, including different loading patterns (flexure or shear), and different material characteristics such as presence of fibers, fibers volume fraction, and presence of shear reinforcement. Moreover, the model is capable of predicting contribution of stirrups and concrete (including compression block and interfacial shear resistance) to shear capacity of UHPFRC beams. • Tensile damage contour predictions along with principal direction, is an effective response parameter for tracing crack propagation zone and failure modes in UHPFRC beams.

•
Removing stirrups does not result in reduction of ductility or load carrying capacity of UHPFRC beams. In other words, UHPFRC beams without shear reinforcement, subjected to dominant shear loading, can attain ultimate moment capacity, without experiencing brittle failure before rebar yielding.