## 1 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 cross-section's flexural capacity. [PERSON] 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. [PERSON] [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. [PERSON] et al. [11] used a cross-section layered model with trilinear tensile stress-strain and stress-crack opening relations to simulate the post-cracking behavior of SFRC. [PERSON] 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. [PERSON] 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, [PERSON] 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.