Flexural Behavior of Steel Fiber Reinforced Concrete Beams: Probabilistic Numerical Modeling and Sensitivity Analysis

: One of the principle issues concerning the practical application of steel ﬁber reinforced concrete (SFRC) is the uncertainty related to its structural behavior, primarily caused by the partially random distribution and orientation of steel ﬁbers in SFRC structural elements. This paper aims to provide a better understanding of how the variance of material properties of the SFRC affects the ﬂexural behavior of SFRC beams. First, a distributed plasticity ﬁber ﬁnite element model of beam ﬂexural behavior is proposed and validated. Then, probability distributions of selected material properties are deﬁned based on existing probabilistic models and experimental results from the literature. Finally, a variance-based sensitivity analysis is performed using Sobol’ indices to identify uncertainties in material properties that contribute most to the uncertainties related to three characteristic points of a beam’s ﬂexural behavior: ﬁrst crack, yield, and collapse point. Sensitivity analysis is performed by surrogating the numerical model using polynomial chaos expansion. The variance in residual tensile strength is identiﬁed as the main contributor to the variance in the ﬂexural behavior of an SFRC beam used in the case study.


Introduction
Brittle failure of concrete in tension can be avoided or delayed by adding randomly dispersed steel fibers to the concrete mixture. Experimental investigations conducted in the last few decades identified several advantages of steel fiber reinforced concrete (SFRC) over conventional reinforced concrete, such as higher tensile and shear strength, improved cyclic response, higher spalling resistance and improved post-cracking behavior [1][2][3][4].
However, SFRC elements exhibit high variance in their structural behavior, mostly due to partially random distribution and orientation of steel fibers. The main factors contributing to this variance can be studied using experimental [5,6] or computational techniques [7,8]. This study focuses on applying computational techniques to investigate how the variance in material properties affects the variance in the flexural behavior of SFRC beams. The basis for such studies are models that can properly simulate the structural behavior of SFRC elements.
Early analytical flexural models of fiber reinforced concrete beams included the effect of steel fibers by modifying concrete's stress-strain relations and then used force equilibrium and the "plane sections remain plane" hypothesis to calculate the crosssection's flexural capacity. Purkiss et al. [9] proposed an analytical model for predicting the load-deformation relation and crack width of beams reinforced with steel bars and steel fibers. To account for the effect of steel fibers, the concrete stress-strain relation in tension is defined as linear with the same elastic modulus as in compression and maximum tensile stress equal to the experimentally obtained modulus of rupture. The stress-strain relation in compression is also modified to include the effect of steel fibers, defined as a function of the fiber aspect ratio, fiber's weight fraction, and experimentally obtained compressive strength. Oh [10] applied the composite materials concept to estimate the flexural strength of SFRC cross-sections with single-and double-sided bar reinforcement. The tensile strength of SFRC is calculated as a sum of the concrete matrix's tensile strength and fibers' strength weighted by their volumetric ratios while considering fibers' orientation, length, and bonding characteristics. Barros et al. [11] used a cross-section layered model with trilinear tensile stress-strain and stress-crack opening relations to simulate the postcracking behavior of SFRC. Qi et al. [12] proposed a flexural strength model for High-Strength-Steel Ultra-High-Performance-Fiber-Reinforced-Concrete (HSS-UHPFRC). Once the maximal tensile strength of UHPFRC is exceeded, the tensile stress drops to the residual tensile strength and stays constant. The residual tensile strength is calculated as the product of the pullout force or fracture failure of a single fiber and the number of fibers crossing a unit crack area, where the fiber failure mode (i.e., pullout or fracture) depends on the fiber aspect ratio. Thereafter, the moment capacity is calculated similarly as in previously mentioned models by imposing an ultimate strain distribution and calculating the moment capacity through cross-section level force equilibrium.
An alternative to analytical models is numerical finite element models that consider the entire structural element. Buzzini et al. [13] simulated the SFRC wall flexural response using a numerical model with solid finite elements and a modified concrete stress-strain relation using the Concrete Damaged Plasticity model as defined in the ABAQUS software, to include the effect of steel fibers on concrete's tensile strength. However, Cunha et al. [7] claim that simulating SFRC as a single-phase material can lead to inaccurate or crude estimates of the behavior of SFRC structural elements, due to different fiber orientation between these elements and test specimens, used to calibrate material behavior laws of SFRC. Therefore, they simulate SFRC as a two-phase material, consisting of a concrete matrix phase and the steel fibers phase. The concrete matrix is simulated using a smeared crack model and solid finite elements, while steel fibers are modeled as embedded elements that can transfer stresses between crack planes of the concrete matrix. A random distribution of fibers inside a structural element is simulated and then inserted into the concrete matrix. Such a model can quantify the variance in SFRC's behavior, using various realizations of potential fiber distributions that lead to different structural behavior that can then be statistically analyzed. Bitencourt Jr. et al. [8] presented a similar numerical model for SFRC composed of three main elements: concrete, concrete-fiber interface, and discrete fibers. The spatial distribution of steel fibers is randomly generated, accounting for the wall effect and meshed independently of concrete. The two are later coupled using a coupling technique for non-aligned meshes. As meso-and micro-scale models are not easily applicable on a structural level, Zhan and Meschke [14] proposed a multilevel modeling framework for SFRC structures. The first level simulates the pullout of a single fiber from the concrete matrix to obtain the pullout force-deformation relation, the second level models the crack bridging effect of a group of fibers using the previously obtained pullout model and the third level represents the structural model that considers the appearance and propagation of cracks and fiber crack-bridging effect using the crack-bridging response obtained from level two. Using such a model, the effect of various design parameters on the structural response of SFRC structures can be investigated starting from the single fiber level.
To summarize, analytical models usually incorporate the uncertainty related to fiber distribution and orientation into several coefficients used to deterministically calculate the tensile and residual tensile strength of SFRC, while more advanced, but computationally more expensive, numerical models explicitly consider this uncertainty using Monte Carlo sampling, to generate an instance of a potential fiber mesh of a structural element. While the randomness related to fiber distribution and orientation might be the dominant source of uncertainty, it is not the only one. Conventional reinforced concrete (i.e., without steel fibers) also has a degree of uncertainty related to its structural behavior as a consequence of the inevitable variance in its material and geometric properties [15]. To the best of the authors' knowledge, none of the currently available studies provide a comparison on the relative contribution of these uncertainties to the uncertainty in the structural behavior of an SFRC element.
This study compares the uncertainties introduced due to steel fibers with the uncertainties inherent in conventional bar reinforced concrete, at different stages of the beam's flexural behavior. The variance of flexural behavior of SFRC beams is examined by conducting a variance-based sensitivity analysis. To that end, a fiber finite element model capable of simulating the flexural behavior of SFRC beams is presented, with the goal of balancing model accuracy and computational efficiency. The proposed model is validated using experimental results available in the literature. Then, uncertainties in the selected material properties are quantified as probability distributions, calibrated according to the available experimental data and existing probabilistic models. Finally, polynomial-chaos-based global sensitivity analysis is used to identify how uncertainties in the material properties contribute to the uncertainties observed at different stages of the flexural behavior of SFRC beams. The results of such variance-based sensitivity analysis can be used to inform variance-reduction actions (i.e., apply computational and/or construction methods focused on decreasing the variance of dominant sources of uncertainty), model dimensionality reduction (i.e., ignore the variance of less influential sources of uncertainty) and can provide qualitative information on how to effectively improve the flexural behavior of SFRC beams by modifying material properties (e.g., how effective is the increase in SFRC's tensile strength in increasing SFRC beam's first crack force).
The structure of the paper is as follows. Section 2 presents methods used for the sensitivity analysis. Section 3 presents the proposed numerical model, adopted stressstrain relations, probability distributions of selected material properties and the definition of first-crack, yield and collapse points. Then, the results of the numerical model validation are presented in Section 4. Section 5 presents a case study of a beam exposed to four-point bending, showing how the proposed numerical model and sensitivity analysis method can be used to identify main sources of uncertainty in different stages of the beam's flexural behavior.

Sensitivity Analysis
A sensitivity analysis of a mathematical model aims to quantify the relative importance of model inputs for a model output. In this paper, the global variance-based sensitivity analysis method (i.e., Sobol' method) [16,17] is used. Model inputs are defined as independent random variables, and Sobol' indices quantify how much the variance of each input or input group contributes to the variance of the model output. It is implied that inputs or input groups that are more important for the model output contribute more to the output's variance and, therefore, have a higher Sobol' index.
The principle idea behind Sobol' indices is to expend the mathematical model into summands of increasing dimensions and then decompose the model output variance into contributions from each summand. If y = f (x) is an integrable function (i.e., a mathematical model), where y is a scalar output and x = x 1 , ..., x n is an input vector of n dimensions, f (x) can be decomposed as follows: where f 0 is a constant and the integrals of the summands with respect to their own variables are zero: Therefore, summands can be calculated by integrating f (x): where x ∼ i stands for vector x that excludes its i-th element. The higher order summands are calculated in the same manner. If f (x) is square-integrable, it follows that all the summands are also square-integrable. Squaring the expansion and integrating over the input domain gives: As the model inputs are random variables, the model output y is also a random variable whose variance (also called total variance) is given by the equation: The total variance can then be decomposed into contributions from each group of inputs (also called partial variances): where following the previous expansion from Equation (6) each partial variance is defined as: Normalizing the partial variances by the total variance, Sobol' indices are obtained: Index S i 1 ,...,i s describes how much of the model output variance is due to the set of model inputs [i 1 , ..., i s ]. The indices are non-negative and sum up to 1. The integer s defines the order of the Sobol' index. For example, S 1 is a first-order Sobol' index showing how much the variance of input 1 alone contributes to the variance of the output. Since one model input can belong to several input groups, total Sobol' indices are used to account for the total effect of an input [18]. The total Sobol' index of an input i represents the sum of all Sobol' indices containing input i: For example, a trivariate function y = f (x 1 , x 2 , x 3 ) can be expended as follows: where can be any univariate, bivariate and trivariate functions of model inputs, respectively, that satisfy Equation (2). The contribution of inputs and input groups to the output variance is quantified using Sobol' indices that sum up to 1: Total Sobol' indices are defined for each input as: The sum of total Sobol' indices of all inputs can be higher than 1, due to the doublecounting of Sobol' indices related to certain input groups.
Monte-Carlo (MC) sampling is usually used to compute Sobol' indices. Saltelli [19] proposed an efficient MC sampling technique to evaluate first-order and total Sobol' indices, which requires a total number of samples of: where N Tot is the total number of samples, N InputDim is the chosen number of samples per model input and n is the number of model inputs. A value often used for N InputDim is 10,000 samples. For the numerical model presented in the next section, which has 6 inputs, the total number of samples would be 80,000. A computationally less demanding alternative to the MC-based Sobol' indices are the polynomial-chaos-expansion-based Sobol' indices [20], which can often be obtained with only a few hundred to a few thousand samples.

Polynomial Chaos Expansion (PCE) Metamodeling
Metamodels (or surrogate models) are, usually, statistically equivalent and computationally less expensive functional approximations of computationally expensive models. The reason for constructing metamodels is to make stochastic modeling (e.g., Monte Carlo simulation) of computationally expensive models feasible. Polynomial Chaos Expansion (PCE) is one of the metamodeling techniques that approximates a computational model as a sum of polynomials orthonormal to the joint probability distribution of the input vector. Advantages of the PCE are that it is a universal metamodel (i.e., does not depend on the type of the computational model), non-intrusive (it is based on repeated runs of the computational model, which is treated as a black box), suited for high-performance computing (i.e., "embarrassingly parallel") and its accuracy converges with the increased number of computational model samples [21].
If y = f (x) is a computational model with an input random vector defined as x = [x 1 , ..., x n ] that follows a joint probabiltiy density function (PDF) p(x), then it's PCE is: where y α are the coefficients (or coordinates) and Φ α (x) are polynomials that are orthonormal with respect to the joint PDF of the input vector p(x). In practice, the infinite sum defined in Equation (18) is truncated to a finite sum, using one of several truncation schemes, such as hyperbolic norm or maximum interaction truncation [22]. If the PCE is truncated to the degree of P, the truncated PCE is: Due to the orthogonality of the polynomial basis with respect to the p(x), the mean and the variance of a PCE are: where A is a subset of coefficients of non-constant polynomial basis elements only.
To estimate the coefficients of the PCE, a computational model is sampled generating a set of input-output pairs, often called experimental design (ED). A commonly used approach for estimating PCE coefficients is the least-squares minimization method [23].

PCE-Based Sobol' Indices
Sudret [20] showed how Sobol' indices of a finite-variance computational model can be obtained by post-processing the coefficients of its PCE metamodel. The principle idea is that the Sobol' decomposition (Equation (1)) of a truncated PCE can be obtained by rearranging the polynomials of a truncated PCE into summands of increasing dimensions: where A i 1 ,...,i s is a set of tuples where only (i 1 , ..., i s ) are non-zero. Comparing the Sobol' decomposition from Equation (1) with Equation (22), the summands of the Sobol' decomposition can be written as: where the polynomial function on the right side depends only on input variables (i 1 , ..., i s ).
Due to the orthonormality of the polynomial basis, the squared integral of the polynomial sum (right side of Equation (22)), as defined in Equation (6), reduces to squaring the coefficients y α . Thus, Sobol' indices of model inputs can be obtained by summing and normalizing squared coefficients of the PCE: Therefore, once the PCE coefficients are obtained, Sobol' indices are calculated analytically with, practically, no computational effort. As obtaining a PCE of a computational model often requires fewer samples than calculating MC-based Sobol' indices, calculating Sobol' indices using the coefficients of a PCE is a viable alternative to the MC-based calculation, assuming that a sufficiently accurate PCE metamodel of the considered computational model can be obtained.

Probabilistic Numerical Model of SFRC Beam Flexural Behavior
Nonlinear flexural behavior of SFRC beams is simulated using a distributed plasticity displacement-based fiber finite element [24] as implemented in OpenSEES [25,26]. Each finite element consists of several cross-sections that serve as integration points, where each cross-section is further discretized into fibers.
To the best of the authors' knowledge, no previous models of SFRC's flexural behavior employed displacement-based fiber finite elements. As this is a distributed plasticity model, the inelastic flexural behavior of the beam is captured along the entire beam length, since each of the fiber cross-sections can simulate the nonlinearities in the SFRC's flexural behavior. Similarly to the analytical models, SFRC is modeled as a single-phase material. However, the uncertainties related to its behavior are considered by defining its material parameters as probability distributions. In terms of the computational cost vs. model accuracy trade-off, the proposed FEM model presents a compromise between simplified computationally inexpensive analytical models that usually consider a single cross-section, and computationally more expensive and complex, but usually more accurate and comprehensive, 3D multi-phase finite element models. Furthermore, fiber-level discretization of cross-sections allows access to stresses and strains in any part of the cross-section along the beam span at any loading level, thus enabling a more detailed insight into the beam behavior without incurring a high computational cost. Due to the small computational cost of such finite element models, probabilistic analysis, such as Monte Carlo sampling and variance-based sensitivity analysis, are feasible.
However, the proposed FEM model has several important limitations. The interaction between shear and flexure is not captured. Shear failure is not considered, hence, only beams that fail in flexure can be simulated using the proposed FEM model. Furthermore, the proposed FEM model may not be able to simulate the deflection-softening behavior of certain beams, as the concrete stress-strain model (Figure 1) does not account for strainsoftening in compression.

Stress-Strain Relations
The effects of steel fibers on concrete mechanical behavior, such as the appearance of residual tensile strength, are accounted for by modifying the stress-strain relation of concrete ( Figure 1). The stress-strain relations used in this study are based on those proposed by [27,28]. Compressive strength, elastic modulus, and the tensile strength of SFRC are assumed to have the same values as for plain concrete. Assuming that the compressive strength of the concrete is known, the equations defined for plain concrete are used to estimate the tensile strength (Equation (25) and the elastic modulus (Equation (28)) [28].
In compression, a bilinear stress-strain relation is used (Figure 1b). Peak compressive strength is estimated to be 85% of the characteristic compressive strength [29].
A trilinear stress-strain relation is used in tension (Figure 1a). Strain corresponding to the peak tensile strength, 1 , can be calculated by dividing the peak tensile strength by the elastic modulus. Residual tensile strength is calculated using Equation (26), assuming that the tensile failure is due to fiber pullout, not fiber fracture. Coefficient a is used to account for the effect that different types of fibers have on the residual tensile strength. For example, due to their geometry, hooked fibers have a higher pullout strength than straight fibers, thus resulting in higher residual tensile strength. The values of the coefficient a are adopted as 1 for straight fibers, 1.5 for irregular fibers, and 2 for hook-end fibers, based on [27]. Strain 2 (Equation (27)) (Figure 1a), after which residual tensile strength is assumed constant, is adopted from [11]. where: V f -volumetric fiber content l d -fiber factor format (i.e., aspect ratio) f ck -characteristic compressive concrete strength 1 -strain at peak tensile strength f t,res -residual tensile strength f t -tensile strength a-coefficient accounting for different types of fiber (a = 1 for straight fibers, a = 1.5 for irregular fibers, a = 2 for hooked fibers) 2 -strain after which residual tensile strength is constant (Figure 1a) Concrete stress-strain relations implemented in OpenSEES could not appropriately represent the proposed SFRC stress-strain relations. Therefore, the Elastic Multi-Linear model, whose nonlinear stress-strain relation is given by a multi-linear curve defined by a set of points, is used [25]. Steel in longitudinal bar reinforcement is modeled using the Giuffré-Menegotto-Pinto model (i.e., Steel02 in OpenSEES). The values of yield strength and elastic modulus varied among experiments and are given in [2,10,12]. The hardening ratio was assumed to be 0.01, and the parameters that control the transition from the elastic to the plastic branch of the steel stress-strain relation were R 0 = 15, cR 1 = 0.925, cR 2 = 0.15 for all specimens used for numerical model validation.

Probability Distribution of Model Parameters
Uncertainties in the material properties are quantified using probability distributions. Table 1 presents the probability distribution types and the assumed Coefficients of Variation (CoV). Mean values of the probability distributions of selected material properties can be obtained experimentally or using Equations (22)- (28). The probability distributions and their CoVs characterizing the uncertainties in concrete compressive and tensile strength, concrete elastic modulus, steel yield strength and steel elastic modulus are adopted based on the recommendations by Wisniweski et al. [30]. It is assumed that the CoVs of material properties defined for plain concrete can be used for SFRC, as the CoVs suggested by [30] are similar or lower to the CoVs of the SFRC's material properties obtained experimentally [2,6].
Residual tensile strength of SFRC depends on the fiber orientation and the number of fibers crossing the crack, both of which are uncertain as they depend on the distribution of fibers along the specimen [31]. Borges et al. [5] and Tiberti et al. [6] experimentally investigated the post-cracking behavior of SFRC and quantified the uncertainties in the residual tensile strength for different crack openings, different displacement application rates, and casting directions. Based on these experimental results, a CoV of 35% is assumed for residual tensile strength. Similarly as for tensile strength, a lognormal distribution is assumed. The variance of geometric properties (e.g., concrete cover or rebar area), steel hardening ratio, fiber content and fiber aspect ratio are not considered in this study. Furthermore, the correlation between considered model inputs is not considered, as the Sobol' method considers models with independent input variables.

Defining Characteristic Points of Beam Flexural Behavior
Three characteristic points in the flexural behavior of beams are identified in this study: the point when the first crack appears (i.e., the first-crack point); the point when longitudinal reinforcement starts to yield (i.e., the yield point) and the point when the beam collapses (i.e., the collapse point). Characteristic points are identified by following the stresses and strains in concrete and longitudinal reinforcement. Figure 2 illustrates how the characteristic points are defined. The appearance of the first crack in the beam is defined as the point when the concrete reaches its residual tensile strength. The yield point is defined as the point when the longitudinal reinforcement reaches its yield strength. Different collapse definitions have been suggested in the literature. Some researchers define collapse by comparing the post-peak and peak loads (e.g., [32,33]). Such a definition is valid for structural elements that exhibit deflection-softening behavior. An alternative is to define maximal allowable strain values. This study uses one such criteria, as proposed by fib Model Code 2010 [34,35]. According to fib Model Code 2020 [34] bending failure (i.e., collapse) of an SFRC beam occurs when one of the following conditions is met: • Attainment of the ultimate compressive strain in the SFRC; • Attainment of the ultimate tensile strain in the longitudinal reinforcement; • Attainment of the ultimate tensile strain in the SFRC.
The ultimate compressive strain for the adopted bilinear concrete stress-strain relation is between 2.4‰and 3.5‰, depending on the concrete compressive strength. The ultimate tensile strain in the longitudinal reinforcement is assumed to be 10% and the ultimate tensile strain in SFRC is equal to 2% for variable strain distribution along the cross-section [34]. Beam shear failure is not considered in this study.

Numerical Model Validation
The proposed numerical model is validated on four-point beam bending experiments available in the literature [2,10,12,36]. All experiments used for validation were performed on a beam reinforced with longitudinal reinforcement and steel fibers, loaded in thirds of the span with two vertical forces of equal intensity. The span-to-depth ratio of the beam was high enough to enable a dominant flexural behavior and all considered beam specimens failed in flexure. Beam parameters that varied among the specimens used for validation were the volumetric ratio of fibers (from 0.38% to 2%), fiber type (straight, irregular and hook-end), fiber aspect ratio (from 30 to 80), concrete strength (from 43 to 135 MPa), longitudinal reinforcement ratio (from 0.67% to 3.4%), beam span (from 1140 mm to 2100 mm) and cross-section dimensions. Details regarding specimens' properties can be found in Table A1 in the Appendix A.
The flexural behavior of the beam is modeled using six finite elements of equal length with two integration points (i.e., cross-sections) per element (Figure 3). Each cross-section is discretized with 20 fibers along the height, representing SFRC and one fiber per bar of longitudinal reinforcement.  Figure 4 compares the load-deflection relation obtained experimentally, using the proposed numerical model for specimen S2V1 tested by Oh [10]. The first crack, yield and collapse points are shown as green, orange and red circles, respectively, showing how the assumed criteria for defining these points compare to the experimentally obtained results. In this case, the numerical model was able to simulate the initial, post-crack, and post-yield stiffness of the beam and to capture the first crack, yield, and collapse points with good accuracy.  Figure 5 compares the experimentally and numerically obtained values of first-crack force for all specimens used for validation. The correlation coefficient between these values is relatively high, at 0.95. Figure 6 presents the relation between the experimentally and numerically obtained values of the first-crack displacement. The correlation coefficient between these values is lower than for the first-crack force, at 0.67. In the case of first-crack displacements, the numerical model almost always underestimated the experimentally obtained displacements. The reason for this underestimation might be that the numerical definition of the first crack point is too conservative since it is defined as the point at which only the bottom fiber of the numerical mid-span cross-section attains residual tensile strength. In the experimental setting, such crack might not have been noticed. For a more significant crack to form, perhaps the residual tensile strength should be attained in several fibers along the height of the cross-section in the numerical model. Nevertheless, the assumed first-crack definition remains on the safe side.  Comparison between numerical and experimental results regarding firstcrack displacement. Figures 7 and 8 show that the numerical model was able to predict the yield force and yield displacement with good accuracy, as the correlation coefficients between the numerical and experimental values are 0.99 and 0.91, respectively.  The correlation coefficient between collapse forces obtained numerically and experimentally is 0.99 (Figure 9). However, the ability of the numerical model to predict collapse displacement is significantly lower, as illustrated in Figure 10. One of the reasons might be that most of the experiments used for validation were load-controlled (red dots in Figure 10), thus the plateau in the load-deflection relation, where the displacement increases significantly, was not captured, underestimating the maximum displacement that the beam could attain. Therefore, load-controlled experiments might not be suitable for validating collapse displacements. Contrary, for all displacement-controlled experiments, the numerical model underestimated the collapse displacement, and thus remained on the safe side (Figure 10), which is most likely the goal of the collapse point criteria used in this paper and defined in fib Model Code 2010 [34].
Furthermore, the identification of the collapse point can be difficult in ductile elements where collapse is not as clearly defined as in brittle elements. The randomness in the structural behavior of SFRC also contributed to the large scatter of experimentally obtained collapse displacements, further reinforcing the need for probabilistic analysis and identification of the most important sources of this uncertainty. Nevertheless, the proposed FEM model can be improved to increase the accuracy of collapse prediction. Potential improvements would be to consider deflection-softening behavior of beams, by, for example, considering strain-softening in compression. Furthermore, other collapse definitions can be considered (e.g., [33]).  Ductility can be expressed by comparing curvatures [32], energies [37] or displacements [38] at ultimate (i.e., collapse) and yield states. In this study, ductility is defined as the ratio between the collapse and the yield displacement, a definition often used for SFRC beams [12,33,39]. Similarly as for collapse displacements, in the case of displacementcontrolled experiments, the numerical model underestimated the ductility of the beams, remaining on the safe side ( Figure 11). Due to the underestimation of beam's collapse displacement in load-controlled experiments, numerical model's estimates of beam ductility often overestimate experimentally obtained values.  Table A2 in the Appendix A provides a summary of first-crack, yield, and collapse forces and displacements obtained experimentally and numerically for all specimens used for validation.

Case Study: SFRC Beam
The RB-SF specimen tested by [2] is used as a case study to illustrate how the proposed numerical model and sensitivity analysis method can be used to identify material properties that are the main sources of uncertainty in the beam's flexural behavior. To illustrate the probabilistic numerical model, 50 realizations of the proposed probabilistic numerical model are performed and the consequent load-deflection relations are presented in Figure 12. The probabilistic model of considered material properties is defined in UQLab software and linked with the numerical model defined in OpenSEESPy using UQLink [26,40,41]. The first-crack, yield, and collapse points are marked with green, orange, and red color, respectively. Higher variance of forces than displacements, related to the three characteristic points, can be observed from the presented results of the probabilistic numerical model, implying that the displacements are less sensitive to the variance of the considered material properties than forces.

PCE Metamodeling
The PCE metamodel of the probabilistic numerical model of the RB-SF specimen used in the previous section is obtained using the UQLab software [41]. Separate PCE metamodels are trained for each of the seven considered Quantities of Interest (QoIs): first-crack displacement and force; yield displacement and force, collapse displacement and force and ductility. An experimental design of 5000 samples is used to train PCEs for all QoIs. Adaptive PCE with a hyperbolic norm truncation scheme is used, and the accuracy of the PCE is measured using the leave-one-out error (LOO error). The details regarding the PCE metamodels for each of the considered QoIs are given in Table 2.

Sensitivity Analysis
Figures 14-20 present the total Sobol' indices for the forces and displacements of three characteristic points of the considered SFRC beams' flexural behavior. Higher total Sobol' index of a material property indicates that the variance of that material property has a higher effect on the variance of the considered QoI, and therefore, changing the value of that material property will have a larger effect on the considered QoI. Material properties with low total Sobol' indices, can be set to a constant value (e.g., to their expected values) without significantly changing the variance of the considered QoI; therefore, the change in the values of material properties with low Sobol' indices will not have a large effect on the considered QoI. This property can be used for dimensionality reduction [42]. Figures 14 and 15 present total Sobol' indices for first-crack force and displacement. The tensile strength of the SFRC has the highest total Sobol' index of around 0.8. The variance of the residual tensile strength contributes around 20% to the variance of firstcrack force and around 10% to the variance of the first-crack displacement. Furthermore, the variance of the elastic modulus of the SFRC is responsible for around 20% of variance of the first-crack displacement but has a very small effect on the first-crack force.    Total Sobol' indices for collapse force and collapse displacement are presented in Figures 18 and 19. Similarly as in the case of yield force, the variance of the collapse force was primarily caused by the variance of the SFRC's residual tensile strength. The variance of the collapse displacement is also dominated by the variance of the SFRC's residual tensile strength. The variance of the SFRC's compressive strength was also important for the variance of the collapse displacement, probably due to certain cases in which beams achieved compression failure, as defined in Section 3.3.   As ductility is defined as the ratio between the collapse and yield displacement, the variance of the ductility depends on the variance of these two displacements. The variance of rebar yield strength and SFRC's residual tensile strength contribute the most to the uncertainty in the SFRC beam's ductility. This means that changing these two material properties is the most effective way to modify the beam's ductility. Thus, the most effective way of increasing the ductility of the SFRC beam considered in this case study is by increasing the residual tensile strength of SFRC (e.g., by increasing the fiber dosage), thereby increasing the collapse displacement by delaying the attainment of 2% strain in tension. However, increasing the residual tensile strength after a certain point will not increase the beam's collapse displacement, as other collapse conditions, such as attainment of maximal strain in compression might be achieved before reaching 2% strain in SFRC tension. Alternatively, the ductility can be increased by decreasing the rebar yield strength, decreasing the beam yield displacement and yield strength, without significantly affecting the beam's collapse displacement ( Figure 18).

Discussion and Conclusions
This paper combines global variance-based sensitivity analysis, (i.e., the Sobol' method) with a fiber finite element model of SFRC beam's flexural behavior to identify material properties whose uncertainty contributes the most to the uncertainty in the flexural behavior of an SFRC beam. To illustrate the proposed approach, a case study SFRC beam is presented. Results of the case study showed that the material parameter whose variance affected the variance of first-crack force and displacement the most (i.e., had the highest Sobol' index) was concrete tensile strength. The variance of the yield point was mostly affected by the variance of residual tensile strength of concrete and yield strength and elastic modulus of the longitudinal reinforcement. Finally, residual tensile strength was the dominant source of uncertainty for the beam's collapse point.
The results of the sensitivity analysis can be used to efficiently reduce the variance of SFRC beams' flexural behavior, by reducing the variance of the material properties with high Sobol' indices. For example, the variance of SFRC's residual tensile strength, primarily caused by the random dispersion of fibers, can be reduced using casting procedures that lead to a more predictive distribution of fibers in the structural element (e.g., horizontal instead of vertical casting), by increasing the fiber content [6] or using advanced computational models to simulate the position and orientation of fibers in an SFRC element, reducing the uncertainty in the distribution of steel fibers along the structural element [43]. Furthermore, the results can help in understanding how changing the values of certain material properties will affect the flexural behavior of the SFRC beam. For example, in the presented case study it was identified that increasing the residual tensile strength will not significantly affect first-crack displacement and force, relevant for serviceability limit-states, but will affect the yield and collapse point.
It should be noted that the results of the Sobol' analysis depend on the assumed probability distributions and their mean and variance values. Therefore, the results presented herein apply for the case of the assumed probability distributions for considered material properties (Table 1). Certain material properties might have low Sobol' indices, due to the assumed low coefficient of variation, not due to their small relevance to the SFRC beam behavior. Furthermore, the limitations of the numerical model, such as its conservative estimation of collapse displacement, should be kept in mind.
Lastly, generalizing the conclusion obtained in this case study to beams of different geometric and material properties requires further analysis. Nevertheless, the presented methodology is applicable to beams of different dimensions and material properties, using the proposed numerical model or more advanced numerical models. Funding: The publication of this research was supported by the company "Construction Center Ltd." from Subotica, Serbia.