Mesoscale Characterization of Fracture Properties of Steel Fiber-Reinforced Concrete Using a Lattice–Particle Model

materials Article Mesoscale Characterization of Fracture Properties of Steel Fiber-Reinforced Concrete Using a Lattice–Particle Model Francisco Montero-Chacón1,*, Héctor Cifuentes 2and Fernando Medina 2 1Department of Engineering, Universidad Loyola Andalucía, Calle Energía Solar 1, 41014 Sevilla, Spain 2ETS de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain; [email protected] (H.C.); [email protected] (F.M.) *Correspondence: [email protected]; Tel.: +34-955-461-600 Academic Editor: Erik Schlangen Received: 16 January 2017; Accepted: 16 February 2017; Published: 21 February 2017 Abstract: This work presents a lattice–particle model for the analysis of steel fiber-reinforced concrete (SFRC). In this approach, fibers are explicitly modeled and connected to the concrete matrix lattice via interface elements. The interface behavior was calibrated by means of pullout tests and a range for the bond properties is proposed. The model was validated with analytical and experimental results under uniaxial tension and compression, demonstrating the ability of the model to correctly describe the effect of fiber volume fraction and distribution on fracture properties of SFRC. The lattice–particle model was integrated into a hierarchical homogenization-based scheme in which macroscopic material parameters are obtained from mesoscale simulations. Moreover, a representative volume element (RVE) analysis was carried out and the results shows that such an RVE does exist in the post-peak regime and until localization takes place. Finally, the multiscale upscaling strategy was successfully validated with three-point bending tests. Keywords: lattice–particle model; fiber-reinforced concrete; fracture 1. Introduction Concrete is today’s most used construction material [ 1 ], and it has been used in many different applications for its versatility and techno-economic advantages. For instance, it allows the construction of structures of any shape and high compressive strength at a fraction of the cost of other important construction materials. On the other hand, its main drawback is its low tensile strength; however, this can be overcome with the inclusion of additional phases such as steel rebars [ 2 ], natural or artificial fibers [3], or carbon nanotubes [4], to cite a few. The demand for high performance concrete (HPC) in civil infrastructures has thus increased and, moreover, empowered the development of ultra-high performance concrete (UHPC) [ 5 ] aiming at extraordinary material features, especially in terms of strength and durability. The inclusion of fibers inside the quasi-brittle concrete matrix enhances the mechanical properties (e.g., strength, toughness, and fatigue life) of the new resulting composite material. In general, the mechanical response of the FRC depends on fiber parameters (size, stiffness, strength, volume content and shape), matrix parameters (stiffness, strength and fracture energy), and fiber–matrix bond parameters (bond strength, stiffness and debonding energy) [ 3 ]. Therefore, it is desirable that the FRC models account for these material parameters. As the use of fiber-reinforced concrete (FRC) becomes more extended, the need for reliable tools to better understand its behavior increases. It is within this context that numerical models play an important role. Certainly, these are less time-demanding (and thus less expensive) than Materials 2017,10, 207; doi:10.3390/ma10020207 www.mdpi.com/journal/materials Materials 2017,10, 207 2 of 19 the experimental campaigns required to characterize FRC mixes. On the other hand, the physical phenomena governing the behavior of FRC must be completely understood so as to solely rely on computational models as a mean for developing new materials in silico, as proposed by the Integrated Computational Materials Engineering (ICME) philosophy [ 6 ]. In any case, experimental tests are an integral part of computational materials science. Certainly, numerical models for the simulation of FRC have experienced important advances. Regarding continuum-based models, different approaches have been followed, such as smeared-crack models with homogenized FRC material properties [ 7 , 8 ] or embedded fibers [ 9 ]; cohesive models [10,11] which can be implemented within an extended-finite element method framework (XFEM); partition of unity finite element method (PUFEM) [ 12 ] which allows for the definition of fiber, matrix, and bond explicit constitutive behaviors; or multiscale approaches such as the micromorphic model [ 13 ] that takes into account the mesostructural level associated to the fiber–matrix interface process at the structural level, or homogenized mesolevel constitutive laws in damage models [14], to cite a few. Discrete models, on the other hand, have shown their suitability for the analysis of fracture mechanics of concrete and other quasi-brittle materials, see for instance refs. [ 15 – 21 ]. For this reason, it has become an appropriate framework for the extension to FRC models, provided fibers can be treated as discrete entities via rigid-body-spring networks [ 22 – 24 ], pure lattice models [ 25 , 26 ] or lattice–particle models [27,28]. In this work, the lattice–particle model developed by the authors [ 28 , 29 ] is enhanced for the analysis of FRC. One main feature of this model is that fibers are explicitly modeled, i.e., the fibers have their own degrees of freedom rather than lumping their arresting effect on the element boundary, as in [ 22 , 27 ]. This is achieved by means of special fiber–matrix interface elements that, in contrast, must be characterized by pullout tests. Moreover, the use of a lattice–particle approach allows for the reduction of the number of degrees of freedom by accounting for different material responses. The model also considers the mechanical and geometrical properties of the constituents of FRC (e.g., mix properties, fiber size and distribution, and material properties) and it has been validated with analytical and experimental results. In this sense, not only the tensile and flexural behaviors have been analyzed, but also the compressive behavior, which has been less studied with numerical models. Finally, the presented model provides a framework for its integration in multiscale analysis of FRC structures. A main concept in multiscale analysis is the existence and definition of a representative volume element (RVE), which has been discussed for plain concrete [ 29 – 31 ]. The effect of fiber reinforcement on the determination of an RVE, especially in the softening regime, is analyzed herein by means of an extensive numerical campaign. 2. Lattice–Particle Model for Fiber-Reinforced Concrete In this section, the main features of the lattice–particle model for the fracture analysis of FRC at the mesoscale (i.e., ~10 mm) are presented. The model is based on previous concepts in [ 17 – 21 ], and it was firstly implemented for plain concrete in [29]. 2.1. Mesostructure Generation At the mesoscale, FRC presents a heterogeneous material structure that includes different phases: mortar, coarse aggregates, interfacial transition zone (ITZ) between these two, and fibers. Therefore, these have to be taken into account in the material generation. In the first place, a coarse aggregate (i.e., maximum aggregate size, d max > 4.75 mm) distribution is generated following the procedure described in [ 20 ], which considers the mix properties (e.g., water to cement (w/c) ratio or aggregates content). A Fuller’s parabola with exponent n= 0.5 is used for the sieve curve [ 32 ] and, as proposed in [ 33 ], a 40% volume content of aggregates is assumed. The generated particles are placed within the domain using the take-and-place method [ 19 , 20 ] and, for the sake of simplicity, spherical shapes are considered. During this process, particle overlaps are not allowed. Materials 2017,10, 207 3 of 19 The main input for the fiber generation is the volume content (V f ), which, in addition to the fiber length (L f ) and diameter (d f ), provides the initial number of one-dimensional inclusions to be placed within the concrete matrix. Fiber elements are introduced at an initial point, within the specimen’s domain, and they are assigned a direction. Both variables are set following a pseudo-random uniform distribution. If the end of the fiber lies out of the boundaries, it is automatically trimmed and the subtracted volume is accounted for when placing the next fiber, preserving the initial V f . Since the aspect ratio of these inclusions is large, no overlap checks are performed at this point. Fibers orientation is biased by the casting direction and rheology of the flow [ 34 ]. Thus, a misorientation angle ( θ ), which controls the maximum angle of the fiber with respect to a specific casting direction, is introduced in order to account for different orientation levels: perfectly oriented fiber distribution with respect to the casting direction ( θ = 0 ◦ ), partially oriented (0 ◦ < θ < 90 ◦ ) or totally misoriented ( θ = 90 ◦ ). Figure 1presents different numerically generated FRC specimens with different Vfand θ. Materials 2017, 10, 207 3 of 18 within the concrete matrix. Fiber elements are introduced at an initial point, within the specimen’s domain, and they are assigned a direction. Both variables are set following a pseudo-random uniform distribution. If the end of the fiber lies out of the boundaries, it is automatically trimmed and the subtracted volume is accounted for when placing the next fiber, preserving the initial Vf. Since the aspect ratio of these inclusions is large, no overlap checks are performed at this point. Fibers orientation is biased by the casting direction and rheology of the flow [34]. Thus, a misorientation angle (θ), which controls the maximum angle of the fiber with respect to a specific casting direction, is introduced in order to account for different orientation levels: perfectly oriented fiber distribution with respect to the casting direction (θ = 0°), partially oriented (0° < θ < 90°) or totally misoriented (θ = 90°). Figure 1 presents different numerically generated FRC specimens with different Vf and θ. (a) (b) (c) Figure 1. FRC numerical mesostrucures: (a) Vf = 0.5% and θ = 90°; (b) Vf = 1.0% and θ = 90°; and (c) Vf = 1.0% and θ = 10°. 2.2. Mesomechanical Elastic Behavior The mechanical model of the current lattice–particle model lies on the interaction between particles and fibers through one-dimensional elements. For this reason, it is necessary to generate a mesh for the concrete matrix, fibers, and fiber–matrix interfaces. The concrete matrix mesh is constructed following a Delaunay’s triangulation using the centroids of the aggregates as in [20] (see Figure 2a). Therefore, the mesh size depends on the aggregates arrangement and every element represents the interaction between two aggregates and their corresponding influence zones, and the area is defined as:   =min , , (1) where ri and rj are the radii of aggregates i and j, respectively. Along with the external nodes of the fibers, some internal nodes are generated in order to link the fibers to the matrix. The fiber element size is defined in terms of the matrix mesh size so as to maintain the average element size of the matrix. In the case of circular cross-section, the area of the fiber elements is defined as:  =  , (2) As mentioned above, the fiber nodes are linked to the closest aggregates through special interface elements, which represent the bond between the concrete matrix and the fibers. The bond elements area is defined as:  = , (3) Every element of the matrix mesh represents the mechanical interaction between the aggregates at the contact point, as depicted in Figure 2b. A set of spring elements acting the in normal and tangential direction, as in [18], is located at this contact point. The normal and tangential stiffness are: Figure 1. FRC numerical mesostrucures: ( a )V f = 0.5% and θ = 90 ◦ ; ( b )V f = 1.0% and θ = 90 ◦ ; and (c)Vf= 1.0% and θ= 10◦. 2.2. Mesomechanical Elastic Behavior The mechanical model of the current lattice–particle model lies on the interaction between particles and fibers through one-dimensional elements. For this reason, it is necessary to generate a mesh for the concrete matrix, fibers, and fiber–matrix interfaces. The concrete matrix mesh is constructed following a Delaunay’s triangulation using the centroids of the aggregates as in [ 20 ] (see Figure 2a). Therefore, the mesh size depends on the aggregates arrangement and every element represents the interaction between two aggregates and their corresponding influence zones, and the area is defined as: Aij =minπr2 i,πr2 j, (1) where riand rjare the radii of aggregates iand j, respectively. Along with the external nodes of the fibers, some internal nodes are generated in order to link the fibers to the matrix. The fiber element size is defined in terms of the matrix mesh size so as to maintain the average element size of the matrix. In the case of circular cross-section, the area of the fiber elements is defined as: Af=πd2 f 4, (2) As mentioned above, the fiber nodes are linked to the closest aggregates through special interface elements, which represent the bond between the concrete matrix and the fibers. The bond elements area is defined as: Ab=πdfLf, (3) Materials 2017,10, 207 4 of 19 Every element of the matrix mesh represents the mechanical interaction between the aggregates at the contact point, as depicted in Figure 2b. A set of spring elements acting the in normal and tangential direction, as in [18], is located at this contact point. The normal and tangential stiffness are: KN ij =α1 Eij Aij Lij , (4) KS ij =α2 Eij Aij Lij , (5) where α1 and α2 are the normal and shear parameters used to adjust the macroscopic elastic modulus and Poisson’s ratio; and their typical values are 1 and 0.2, respectively. On the other hand, E ij is the local (i.e., element) elastic modulus and L ij is the length of the element directly obtained from the triangulation. The local elastic modulus is defined in terms of the three phases that are considered in the concrete matrix, mortar (E m ), coarse aggregates (E a ), and ITZ (E ITZ ), as well as on the aggregate radii, the ITZ thickness of each aggregate (L iITZ and L jITZ ), and the mortar length (L ijm ), which is the rest of the element. Assuming a serial coupling along the element [12], the local elastic modulus reads: Lij Eij =ri Ea +LITZ i EITZ +Lm ij Em +LITZ j EITZ +rj Ea, (6) The elastic modulus for cement mortars varies from 10 GPa to 70 GPa, depending on the porosity [ 35 ]; and that of coarse aggregates ranges from 70 GPa to 90 GPa [ 35 , 36 ]. For the numerical simulations, the material properties considered are obtained from [ 20 ], although lower scale discrete models can be also used to compute these [14]. Fibers are modeled as truss elements, with the following axial stiffness: KN f=EfAf Lf,e , (7) where E f is the elastic modulus of the fiber, and A f the cross-sectional area of the fiber, given by Equation (2). The bond elements (Figure 2c) are defined as pure shear elements, since no axial interaction is expected at the interface: KS b=EbAb Lb , (8) In general, the normal and tangential displacements of an element, un and us , respectively, can be computed as: un=(u·n)n, (9) us=u−un, (10) where u is the element displacement vector and n is the normal vector that, for an element connecting nodes iand j, is defined as: n=xj−xi Lij , (11) with xiand xjthe coordinates of nodes iand j, respectively. The normal and tangential element strains can be obtained as: εn=|un| Lij , (12) Materials 2017,10, 207 5 of 19 εs=|us| Lij , (13) and, finally, the normal and shear stresses can respectively be computed from: σ=α1Eεn, (14) τ=α2Eεs, (15) Materials 2017, 10, 207 5 of 18 (a) (b) (c) Figure 2. Mesomechanical features of the lattice–particle model: (a) concrete matrix Delaunay’s triangulation and fiber–matrix interaction; (b) aggregates interaction at the contact point and spring elements; and (c) resulting matrix, bond, and fiber mesh. 2.3. Mesoscale Fracture Behavior In order to account for the material non-linearity intrinsic to fracture, a sequentially event-driven solution scheme [37], in which the non-linear problem is solved stepwise into several linear analyses, is used. Thus, at every simulation step, unitary prescribed forces (or displacements) are applied to the specimen and the element stress pairs are calculated. In the case of the concrete matrix, a linear softening curve is considered (although other shapes can be implemented). This law is discretized into N segments (see Figure 3a), and the elastic modulus of the element with the maximum stress-to-strength ratio is updated following this expression: = ,   󰇛󰇜   , (16) where i is the current segment; E i and f t,i are the current elastic modulus and tensile strength, respectively; ε t is the initial cracking strain; and ε f is the ultimate tensile strain which is defined in terms of the element length and the fracture energy, G F , to avoid the mesh-dependency issues: =    +    , (17) As observed in Figure 3a, with the discretization of the softening curve, some fracture energy is lost. Such energy loss can be minimized by increasing the number of segments, N. On the other hand, a corrective term  is included such that an updated value for the fracture energy is used, ′=, which is equivalent to increasing the tensile strength and ultimate strain  times, as proposed in [37]. As a reference, for a total number of segments N = 10, the corrective term yields ~1.1. In the case of bond and fiber elements, a general bilinear material behavior is implemented so as to account for a wide range of behaviors, namely from softening to hardening. This is achieved by defining initial and ultimate stress–strain pairs, and the sequential reduction takes place therein, as explained above. Since shear interaction is present in matrix and bond elements, a Mohr–Coulomb failure surface with tension cut-off and compression cap (Figure 3b) is used, as proposed in [17]. This surface is defined by the tensile strength, f t ; the slope of the failure envelope or angle of internal friction, ϕ; the cohesion stress, c; and the compressive strength, f c . Although the ITZ is known to be the weakest link in the concrete mesostructure, it has little influence on the macroscopic tensile strength [38], in contrast to mortar properties. Therefore, local fracture parameters in the matrix elements are those for mortar. In this regard, tensile strength in common mortars vary from 2 MPa up to 15 MPa [36]. For granular materials, internal friction angles are in the range of 25° to 45° [39]. The cohesion stress is set as c = 2f t . It must be remarked that, Figure 2. Mesomechanical features of the lattice–particle model: ( a ) concrete matrix Delaunay’s triangulation and fiber–matrix interaction; ( b ) aggregates interaction at the contact point and spring elements; and (c) resulting matrix, bond, and fiber mesh. 2.3. Mesoscale Fracture Behavior In order to account for the material non-linearity intrinsic to fracture, a sequentially event-driven solution scheme [37], in which the non-linear problem is solved stepwise into several linear analyses, is used. Thus, at every simulation step, unitary prescribed forces (or displacements) are applied to the specimen and the element stress pairs are calculated. In the case of the concrete matrix, a linear softening curve is considered (although other shapes can be implemented). This law is discretized into Nsegments (see Figure 3a), and the elastic modulus of the element with the maximum stress-to-strength ratio is updated following this expression: Ei=ft,i εt−(i−1)εf−εt N, (16) where iis the current segment; E i and f t,i are the current elastic modulus and tensile strength, respectively; εt is the initial cracking strain; and εf is the ultimate tensile strain which is defined in terms of the element length and the fracture energy, GF, to avoid the mesh-dependency issues: εf=2GF ftL+ft E0, (17) As observed in Figure 3a, with the discretization of the softening curve, some fracture energy is lost. Such energy loss can be minimized by increasing the number of segments, N. On the other hand, a corrective term ς is included such that an updated value for the fracture energy is used, G0F=ς2GF , which is equivalent to increasing the tensile strength and ultimate strain ς times, as proposed in [ 37 ]. As a reference, for a total number of segments N= 10, the corrective term yields ς∼1.1. In the case of bond and fiber elements, a general bilinear material behavior is implemented so as to account for a wide range of behaviors, namely from softening to hardening. This is achieved by defining initial and ultimate stress–strain pairs, and the sequential reduction takes place therein, as explained above. Materials 2017,10, 207 6 of 19 Since shear interaction is present in matrix and bond elements, a Mohr–Coulomb failure surface with tension cut-off and compression cap (Figure 3b) is used, as proposed in [ 17 ]. This surface is defined by the tensile strength, f t ; the slope of the failure envelope or angle of internal friction, ϕ ; the cohesion stress, c; and the compressive strength, fc. Although the ITZ is known to be the weakest link in the concrete mesostructure, it has little influence on the macroscopic tensile strength [ 38 ], in contrast to mortar properties. Therefore, local fracture parameters in the matrix elements are those for mortar. In this regard, tensile strength in common mortars vary from 2 MPa up to 15 MPa [ 36 ]. For granular materials, internal friction angles are in the range of 25 ◦ to 45 ◦ [ 39 ]. The cohesion stress is set as c= 2f t . It must be remarked that, although this mortar information is lumped into the spring sets, it could be explicitly modeled by including mortar particles, and thus decreasing the dmax value of the generated distribution. Materials 2017, 10, 207 6 of 18 although this mortar information is lumped into the spring sets, it could be explicitly modeled by including mortar particles, and thus decreasing the d max value of the generated distribution. (a) (b) Figure 3. Mesoscale fracture behavior of matrix elements: (a) linear softening curve with sequential reduction; and (b) Mohr–Coulomb fracture surface with tension cut-off and compression cap. Even though fracture events are more likely to happen in matrix and bond elements, due to their lower strength, fibers are assigned a Rankine failure criterion characterized by the fiber yield strength. 2.4. Meso-Macro Upscaling Strategy The fiber-reinforced lattice–particle model can be used to upscale the mesoscale properties to the macroscale (i.e., ~1 m), in a numerical homogenization fashion [40]. However, in the case of quasi-brittle materials, that show strain localization, this step is not straightforward and some considerations must be taken into account. One main concept in homogenization-based schemes is the so-called representative volume element (RVE) that, basically, is a volume element over which a measured value becomes representative of the whole. In plain concrete, the definition of such an RVE is not possible due to strain localization [29,30], but can be overcome by the definition of failure average zones [31]. However, FRC may present in some cases a pseudo-plastic behavior, and thus circumventing this problem. It is within this context that numerical homogenization can be used as a multiscale approach so as to evaluate the macroscopic behavior of FRC structures. In this work we propose the discrete to continuum upscaling procedure presented in [14], by which the material properties of the continuum damage plasticity (CDP) model are obtained by virtual testing. Thus, the tensile cracking strain, which is required by the CDP model to describe the tensile behavior of the material, can be obtained as: 0   ck t tt , E   (18) where ε t is the total tensile strain, σ t is the tensile stress and E 0 is the elastic modulus, obtained by means of mesoscale simulations. 3. Results and Discussion The model presented in Section 2 considers the matrix, fiber, and fiber–matrix interaction parameters required to describe the mechanical behavior of FRC, as suggested in [3]. Therefore, it requires the input of 15 material properties and the material specifications (e.g., mix information, and fiber distribution). These properties are summarized in Table 1. In this section, the fiber-reinforced lattice particle model is calibrated and validated by means of different tests. In the first place, pullout tests are carried out to characterize the fiber–matrix interface behavior, which is essential in the response of FRC models. Next, uniaxial tensile and compressive tests are carried out so as to validate the model. Finally, the upscaling approach described in Section Figure 3. Mesoscale fracture behavior of matrix elements: ( a ) linear softening curve with sequential reduction; and (b) Mohr–Coulomb fracture surface with tension cut-off and compression cap. Even though fracture events are more likely to happen in matrix and bond elements, due to their lower strength, fibers are assigned a Rankine failure criterion characterized by the fiber yield strength. 2.4. Meso-Macro Upscaling Strategy The fiber-reinforced lattice–particle model can be used to upscale the mesoscale properties to the macroscale (i.e., ~1 m), in a numerical homogenization fashion [ 40 ]. However, in the case of quasi-brittle materials, that show strain localization, this step is not straightforward and some considerations must be taken into account. One main concept in homogenization-based schemes is the so-called representative volume element (RVE) that, basically, is a volume element over which a measured value becomes representative of the whole. In plain concrete, the definition of such an RVE is not possible due to strain localization [29,30] , but can be overcome by the definition of failure average zones [ 31 ]. However, FRC may present in some cases a pseudo-plastic behavior, and thus circumventing this problem. It is within this context that numerical homogenization can be used as a multiscale approach so as to evaluate the macroscopic behavior of FRC structures. In this work we propose the discrete to continuum upscaling procedure presented in [ 14 ], by which the material properties of the continuum damage plasticity (CDP) model are obtained by virtual testing. Thus, the tensile cracking strain, which is required by the CDP model to describe the tensile behavior of the material, can be obtained as: e εck t=εt−σt E0, (18) where εt is the total tensile strain, σt is the tensile stress and E 0 is the elastic modulus, obtained by means of mesoscale simulations. Materials 2017,10, 207 7 of 19 3. Results and Discussion The model presented in Section 2considers the matrix, fiber, and fiber–matrix interaction parameters required to describe the mechanical behavior of FRC, as suggested in [ 3 ]. Therefore, it requires the input of 15 material properties and the material specifications (e.g., mix information, and fiber distribution). These properties are summarized in Table 1. In this section, the fiber-reinforced lattice particle model is calibrated and validated by means of different tests. In the first place, pullout tests are carried out to characterize the fiber–matrix interface behavior, which is essential in the response of FRC models. Next, uniaxial tensile and compressive tests are carried out so as to validate the model. Finally, the upscaling approach described in Section 2.4 is validated via three-point bending tests (3PBT), showing the ability of the model to reproduce the fracture behavior FRC. Table 1. Summary of the input parameters required by the fiber-reinforced lattice–particle model. Phase Input Properties Matrix Elastic: Em,Ea Fracture: ft,fc,c,ϕ,GF Mix information: w/c,a/c,d max Fiber Elastic: Ef Fracture: σ1,σ2,εf/εi Bond Elastic: Eb Fracture: τ1,τ2,εf/εi 3.1. Pullout Test The fiber–matrix interface behavior is best characterized by means of fiber pullout tests [ 3 ]. In these tests, the fiber is partially embedded in the matrix and it is pulled out through its free end, as depicted in Figure 4. The global response is governed by the so-called force-slip curve, and three main stages can be observed in this curve: (1) an elastic branch prior to the chemical strength degradation; (2) non-linear fiber–matrix debonding; and (3) frictional softening [3]. Materials2017,10,207 7of18 2.4isvalidatedviathree‐pointbendingtests(3PBT),showingtheabilityofthemodeltoreproduce thefracturebehaviorFRC. Table1.Summaryoftheinputparametersrequiredbythefiber‐reinforcedlattice–particlemodel. PhaseInputProperties Matrix Elastic:Em,Ea Fracture:ft,fc,c,φ,GF Mixinformation:w/c,a/c,dmax FiberElastic:Ef Fracture:σ1,σ2,εf/εi BondElastic:Eb Fracture:τ1,τ2,εf/εi 3.1.PulloutTest Thefiber–matrixinterfacebehaviorisbestcharacterizedbymeansoffiberpullouttests[3].In thesetests,thefiberispartiallyembeddedinthematrixanditispulledoutthroughitsfreeend,as depictedinFigure4.Theglobalresponseisgovernedbytheso‐calledforce‐slipcurve,andthree mainstagescanbeobservedinthiscurve:(1)anelasticbranchpriortothechemicalstrength degradation;(2)non‐linearfiber–matrixdebonding;and(3)frictionalsoftening[3].  (a)(b) Figure4.Fiberpullouttestnumericalsimulationsetup:(a)thefiber(incolorblue)ispartially embeddedintheconcretematrixanditispulledoutbyitsfreeend;and(b)thepullouttestallows forthecharacterizationoftheinterfaceelements(incolorred). Thepresentmodelaccountsfordifferentforce‐slipbehaviors:sliphardening,perfectly elastic‐plastic,andslipsoftening.Theseareimplementedviabilinearlawsthatrequirefourmaterial properties,asdescribedinTable1. SeveralpullouttestsweresimulatedandcomparedtotheresultsbyCunhaetal.[41]and Kimetal.[42].Inthesesteelfiberpulloutsimulations,differentconcretespecimenswithdifferent fibersizeswereconsidered,inaccordancewith[41,42].Thespecimenswerefixedatthebottom surfaceandthefreeendofthehalf‐embeddedsteelfiberwaspulled‐out. Figure5showstheresultsforthepulloutsimulations.Ascanbeobserved,themodelisableto reproducewellthepulloutforcesatlowslipvalues,i.e.,below1mm(Figure5a).Ontheotherhand, themodelisalsoabletodeveloplargeslipvalues(Figure5b).Table2presentstheoptimalrangesfor theinterfacematerialpropertiesthatareusedinfurthersimulations. Figure 4. Fiber pullout test numerical simulation setup: ( a ) the fiber (in color blue) is partially embedded in the concrete matrix and it is pulled out by its free end; and ( b ) the pullout test allows for the characterization of the interface elements (in color red). The present model accounts for different force-slip behaviors: slip hardening, perfectly elastic-plastic, and slip softening. These are implemented via bilinear laws that require four material properties, as described in Table 1. Several pullout tests were simulated and compared to the results by Cunha et al. [ 41 ] and Kim et al. [ 42 ]. In these steel fiber pullout simulations, different concrete specimens with different fiber Materials 2017,10, 207 8 of 19 sizes were considered, in accordance with [ 41 , 42 ]. The specimens were fixed at the bottom surface and the free end of the half-embedded steel fiber was pulled-out. Figure 5shows the results for the pullout simulations. As can be observed, the model is able to reproduce well the pullout forces at low slip values, i.e., below 1 mm (Figure 5a). On the other hand, the model is also able to develop large slip values (Figure 5b). Table 2presents the optimal ranges for the interface material properties that are used in further simulations. Materials 2017, 10, 207 8 of 18 (a) (b) Figure 5. Fiber pullout test numerical simulation results: (a) low slip range; and (b) high slip range. Table 2. Fiber–matrix interface material properties range. Property Values Eb (GPa) 1–10 τ1 (MPa) 1–10 τ2 (MPa) <2 εf/εi 5–15 3.2. Tensile Tests Uniaxial tensile tests provide important information regarding the fracture properties of FRC, for instance the tensile strength (f t ) and fracture energy (G F ). Moreover, the effect of fiber volume content and distribution can be observed. In this work, the uniaxial tensile response of steel fiber-reinforced concrete (SFRC) is validated using the analytical model proposed by Karihaloo and Lange-Kornbak [43] and experimental double-edge-notched (DEN) tension tests by Gopalaratnam and Shah [44]. On the other hand, the influence of V f and θ in the tensile response is analyzed by means of unnotched uniaxial tension tests. 3.2.1. Elastic Modulus In the first place, and before characterizing the fracture behavior of SFRC, the elastic modulus of the composite material was determined following an automatic numerical campaign, evaluating a total of 100 prismatic specimens with size 50 × 25 × 150 mm 3 , d max = 8 mm, steel fibers with L f = 35 mm and d f = 0.5 mm, random distribution and variable V f = 0.0%–3.0%. The elastic properties of the phases were E m = 20 GPa, E a = 70 GPa, E f = 210 GPa, and E b = 5 GPa. In Figure 6, the numerical results are compared to the analytical expression from [43] and the rule of mixture bounds for composite materials. It can be observed that the numerical results are in agreement with those obtained from [43] and the trend is similar to that of the parallel behavior, as Figure 5. Fiber pullout test numerical simulation results: (a) low slip range; and (b) high slip range. Table 2. Fiber–matrix interface material properties range. Property Values Eb(GPa) 1–10 τ1(MPa) 1–10 τ2(MPa) <2 εf/εi5–15 3.2. Tensile Tests Uniaxial tensile tests provide important information regarding the fracture properties of FRC, for instance the tensile strength (f t ) and fracture energy (G F ). Moreover, the effect of fiber volume content and distribution can be observed. In this work, the uniaxial tensile response of steel fiber-reinforced concrete (SFRC) is validated using the analytical model proposed by Karihaloo and Lange-Kornbak [ 43 ] and experimental double-edge-notched (DEN) tension tests by Gopalaratnam and Shah [ 44 ]. On the other hand, the influence of V f and θ in the tensile response is analyzed by means of unnotched uniaxial tension tests. Materials 2017,10, 207 9 of 19 3.2.1. Elastic Modulus In the first place, and before characterizing the fracture behavior of SFRC, the elastic modulus of the composite material was determined following an automatic numerical campaign, evaluating a total of 100 prismatic specimens with size 50 × 25 × 150 mm 3 ,d max = 8 mm, steel fibers with Lf= 35 mm and d f = 0.5 mm, random distribution and variable V f = 0.0%–3.0%. The elastic properties of the phases were Em= 20 GPa, Ea= 70 GPa, Ef= 210 GPa, and Eb= 5 GPa. In Figure 6, the numerical results are compared to the analytical expression from [ 43 ] and the rule of mixture bounds for composite materials. It can be observed that the numerical results are in agreement with those obtained from [ 43 ] and the trend is similar to that of the parallel behavior, as expected. On the other hand, it can be also observed that as V f increases, the scatter decreases as a result of an increase in the number of nodes in the system. Materials 2017, 10, 207 9 of 18 expected. On the other hand, it can be also observed that as V f increases, the scatter decreases as a result of an increase in the number of nodes in the system. Figure 6. Effect of V f in the elastic modulus of SFRC. 3.2.2. Tensile Strength A total number of 50 specimens with the same dimensions and properties used for the elastic modulus analysis, were subjected to complete uniaxial tension tests. The fracture properties for the concrete matrix were f t = 4 MPa, f c = −12f t , c = 2f t , ϕ = 35°, and G F = 10 N/m; and a softening behavior for the bond determined by τ 1 = 1 MPa, τ 2 = 0.1 MPa, and ε f /ε i = 15 was considered. The volume fraction range was set to V f = 0.5%–1.5% so as to analyze the experimental configurations in more details. Figure 7 shows the numerical results for the tensile strength for different volume fractions. The results are in agreement with those obtained from the expression in [43]. The effect of V f in f t is slightly lower in the numerical simulations, but larger than those reported in [45]. Figure 7. Effect of V f in the tensile strength of SFRC. Splitting tensile tests were carried out on cylindrical specimens of 100 mm in the laboratory, using the same concrete mix and fiber geometry as the simulations, and a value of f st = 4.97 MPa was obtained for V f = 0.6%. The equivalent tensile strength, obtained with the Model Code formula [46], is thus f t = 4.47 MPa, which is in agreement with the numerical result, f t = 4.21 MPa. The previous simulations serve as the basis for the analysis of the effect of V f and fiber orientation in the stress–strain uniaxial tension response. Figure 8 presents the effect on the tensile response of V f in random distribution (θ = 60°) and fiber orientation (for a given V f = 1.0%). In the first case (Figure 8a), as V f increases, the tensile strength and ductility increases. In the second case (Figure 8b), as the fibers become more aligned with respect to the loading axis, the strength and ductility increases. It must be remarked that this effect is more evident for high θ values; however, the tortuosity of the crack may diminish this effect for highly oriented distributions. In Figure 8b, it can be observed that for misorientations below 15°, the maximum arresting effect of the fibers is approached. Figure 6. Effect of Vfin the elastic modulus of SFRC. 3.2.2. Tensile Strength A total number of 50 specimens with the same dimensions and properties used for the elastic modulus analysis, were subjected to complete uniaxial tension tests. The fracture properties for the concrete matrix were f t = 4 MPa, f c = − 12f t ,c= 2f t , ϕ = 35 ◦ , and G F = 10 N/m; and a softening behavior for the bond determined by τ1 = 1 MPa, τ2 = 0.1 MPa, and εf / εi = 15 was considered. The volume fraction range was set to V f = 0.5%–1.5% so as to analyze the experimental configurations in more details. Figure 7shows the numerical results for the tensile strength for different volume fractions. The results are in agreement with those obtained from the expression in [ 43 ]. The effect of V f in f t is slightly lower in the numerical simulations, but larger than those reported in [45]. Materials 2017, 10, 207 9 of 18 expected. On the other hand, it can be also observed that as V f increases, the scatter decreases as a result of an increase in the number of nodes in the system. Figure 6. Effect of V f in the elastic modulus of SFRC. 3.2.2. Tensile Strength A total number of 50 specimens with the same dimensions and properties used for the elastic modulus analysis, were subjected to complete uniaxial tension tests. The fracture properties for the concrete matrix were f t = 4 MPa, f c = −12f t , c = 2f t , ϕ = 35°, and G F = 10 N/m; and a softening behavior for the bond determined by τ 1 = 1 MPa, τ 2 = 0.1 MPa, and ε f /ε i = 15 was considered. The volume fraction range was set to V f = 0.5%–1.5% so as to analyze the experimental configurations in more details. Figure 7 shows the numerical results for the tensile strength for different volume fractions. The results are in agreement with those obtained from the expression in [43]. The effect of V f in f t is slightly lower in the numerical simulations, but larger than those reported in [45]. Figure 7. Effect of V f in the tensile strength of SFRC. Splitting tensile tests were carried out on cylindrical specimens of 100 mm in the laboratory, using the same concrete mix and fiber geometry as the simulations, and a value of f st = 4.97 MPa was obtained for V f = 0.6%. The equivalent tensile strength, obtained with the Model Code formula [46], is thus f t = 4.47 MPa, which is in agreement with the numerical result, f t = 4.21 MPa. The previous simulations serve as the basis for the analysis of the effect of V f and fiber orientation in the stress–strain uniaxial tension response. Figure 8 presents the effect on the tensile response of V f in random distribution (θ = 60°) and fiber orientation (for a given V f = 1.0%). In the first case (Figure 8a), as V f increases, the tensile strength and ductility increases. In the second case (Figure 8b), as the fibers become more aligned with respect to the loading axis, the strength and ductility increases. It must be remarked that this effect is more evident for high θ values; however, the tortuosity of the crack may diminish this effect for highly oriented distributions. In Figure 8b, it can be observed that for misorientations below 15°, the maximum arresting effect of the fibers is approached. Figure 7. Effect of Vfin the tensile strength of SFRC. Materials 2017,10, 207 16 of 19 Materials 2017, 10, 207 15 of 18 3.4.2. Experimental Validation Three-point bending tests (3PBT) were carried out on prismatic beams with a cross section of 100 mm × 100 mm, 400 mm of span, and a total length of 440 mm. The FRC composition was the same as in Section 3.3. The beams were notched with a thin (3 mm) diamond saw to a notch to depth ratio a/W = 1/6, according to EN 14651. The crack mouth opening displacement (CMOD) was measured with a clip gauge transducer and used as the feedback control signal. The load-point deflection was measured simultaneously by means of a linearly variable displacement transducer (LVDT) mounted on a rigid frame in order to avoid parasitic torsional effects on the measurement of vertical displacements. The tests were performed in a stiff closed-loop universal testing machine with a maximum load capacity of 50 kN. The 3PBT numerical simulations were carried out in Abaqus using the concrete damage plasticity model [48], which has been successfully proven in these kind of analyses [14]. Figure 15a shows the experimental setup for the 3PBT. In 3PBT configurations, tensile behavior governs the global response of the structure, therefore the uniaxial tensile behavior of the material was obtained from lattice–particle simulations on specimens fulfilling the previous RVE size considerations, and the resulting stress–cracking strain and tensile damage–cracking strain curve were extracted (Figure 15b). (a) (b) Figure 15. (a) 3PBT experimental setup; and (b) uniaxial tensile stress–cracking strain and tensile damage–cracking strain curve obtained for V f = 0.6% and θ = 30°. Figure 16 presents the load-CMOD response comparison between the experimental and numerical results. The numerical results are in agreement with the experimental results, especially in the peak and post-peak range, showing the ability of the mesoscale fiber-reinforced lattice–particle model for providing reliable input parameters for macroscale models in a hierarchical homogenization multiscale approach. Figure 16. 3PBT load-CMOD validation, V f = 0.6% and θ = 30°. Figure 15. ( a ) 3PBT experimental setup; and ( b ) uniaxial tensile stress–cracking strain and tensile damage–cracking strain curve obtained for Vf= 0.6% and θ= 30◦. Figure 16 presents the load-CMOD response comparison between the experimental and numerical results. The numerical results are in agreement with the experimental results, especially in the peak and post-peak range, showing the ability of the mesoscale fiber-reinforced lattice–particle model for providing reliable input parameters for macroscale models in a hierarchical homogenization multiscale approach. Materials 2017, 10, 207 15 of 18 3.4.2. Experimental Validation Three-point bending tests (3PBT) were carried out on prismatic beams with a cross section of 100 mm × 100 mm, 400 mm of span, and a total length of 440 mm. The FRC composition was the same as in Section 3.3. The beams were notched with a thin (3 mm) diamond saw to a notch to depth ratio a/W = 1/6, according to EN 14651. The crack mouth opening displacement (CMOD) was measured with a clip gauge transducer and used as the feedback control signal. The load-point deflection was measured simultaneously by means of a linearly variable displacement transducer (LVDT) mounted on a rigid frame in order to avoid parasitic torsional effects on the measurement of vertical displacements. The tests were performed in a stiff closed-loop universal testing machine with a maximum load capacity of 50 kN. The 3PBT numerical simulations were carried out in Abaqus using the concrete damage plasticity model [48], which has been successfully proven in these kind of analyses [14]. Figure 15a shows the experimental setup for the 3PBT. In 3PBT configurations, tensile behavior governs the global response of the structure, therefore the uniaxial tensile behavior of the material was obtained from lattice–particle simulations on specimens fulfilling the previous RVE size considerations, and the resulting stress–cracking strain and tensile damage–cracking strain curve were extracted (Figure 15b). (a) (b) Figure 15. (a) 3PBT experimental setup; and (b) uniaxial tensile stress–cracking strain and tensile damage–cracking strain curve obtained for V f = 0.6% and θ = 30°. Figure 16 presents the load-CMOD response comparison between the experimental and numerical results. The numerical results are in agreement with the experimental results, especially in the peak and post-peak range, showing the ability of the mesoscale fiber-reinforced lattice–particle model for providing reliable input parameters for macroscale models in a hierarchical homogenization multiscale approach. Figure 16. 3PBT load-CMOD validation, V f = 0.6% and θ = 30°. Figure 16. 3PBT load-CMOD validation, Vf= 0.6% and θ= 30◦. 4. Conclusions In this work, a lattice–particle model for the analysis of fracture properties of steel fiber-reinforced concrete has been presented. The model has been validated with respect to existing analytical models [ 43 ] and experimental data [ 44 ], in the uniaxial tensile behavior. Moreover, an experimental campaign was carried out in order to validate the compressive and flexural behaviors. The characterization of the fiber–matrix interface is of great importance in the presented model. For this reason, pullout tests were carried out and compared to existing experimental results [ 41 , 42 ], and the ability of the model to account for different bonding behaviors has been shown. In this sense, a range for the bonding material properties has been presented. The fiber-reinforced lattice–particle model is able to account for different volume fractions (V f ) and fiber orientations, characterized by the fiber misorientation angle ( θ ), and the results showed that an increase in V f leads to an increase in the ductility, while an increase in θ has the opposite effect. The model was used to analyze the effect of such parameters on the fracture properties and Materials 2017,10, 207 17 of 19 the results agreed well with the analytical and experimental results, especially for the tensile and compressive strengths. The effect on the ductility was also analyzed by means of the fracture energy and, although the trend is similar to the analytical model by Karihaloo and Lange-Kornbak [ 43 ], the predicted values were lower. An upscaling hierarchical homogenization-based scheme has been presented in order to provide material information at the macroscale by performing virtual tests on the mesoscale model. In this regard, a representative volume element analysis was carried out for FRC in order to determine its size. Moreover, an RVE can be found for the hardening and softening regimes, at least until localization takes place. Thus, conventional multiscale approaches may be followed without loss of generality if these conditions are satisfied. The ability of the lattice–particle model to provide material properties within a multiscale framework has been thus demonstrated by means of 3PBT. Finally, not only can the presented model be used to provide material properties at the mesolevel, but it can also be connected to other physical models which provide material properties at lower (or larger scales) or information of the material structure, e.g., casting process modeling, opening the door to an Integrated Computational Materials Engineering (ICME) approach [ 6 ] for the design of fiber-reinforced cement-based composites. Acknowledgments: The authors would like to acknowledge the financial support from the research project BIA2013-48352-P (Ministry of Economy and Competitiveness of Spain). Author Contributions: Francisco Montero-Chacón and Fernando Medina conceived the lattice–particle model and performed the numerical simulations and Héctor Cifuentes performed the experimental campaign and all the authors wrote the article. Conflicts of Interest: The authors declare no conflict of interest. References 1. Kosmatka, S.; Panarese, W.; Kerkhoff, B. Design and Control of Concrete Mixtures; Portland Cement Association: New York, NY, USA, 2011. 2. MacGregor, J.; Wight, J.; Teng, S.; Irawan, P. Reinforced Concrete: Mechanics and Design; Prentice Hall: Upper Saddle River, NJ, USA, 1997. 3. Bentur, A.; Mindess, S. Fibre Reinforced Cementitious Composites; Taylor & Francis: London, UK, 2007. 4. Konsta-Gdoutos, M.; Metaxa, Z.; Shah, S. Highly dispersed carbon nanotube reinforced cement based materials. Cem. Concr. Res. 2010,40, 1052–1059. [CrossRef] 5. Wang, C.; Yang, C.; Liu, F.; Wan, C.; Pu, X. Preparation of ultra-high performance concrete with common technology and materials. Cem. Concr. Comp. 2012,34, 538–544. [CrossRef] 6. Horstemeyer, M. Integrated Computational Materials Engineering (ICME) for Metals: Using Multiscale Modeling to Invigorate Engineering Design with Science; John Wiley & Sons: Hoboken, NJ, USA, 2012. 7. Barros, J.A.; Figueiras, J.A. Model for the analysis of steel fibre reinforced concrete slabs on grade. Comput. Struct. 2001,79, 97–106. [CrossRef] 8. Özcan, D.M.; Bayraktar, A.; Sahin, A.; Haktanir, T.; Türker, T. Experimental and finite element analysis on the steel fiber-reinforced concrete (SFRC) beams ultimate behavior. Constr. Build. Mater. 2009 ,23, 1064–1077. [CrossRef] 9. Cunha, V.M.C.F.; Barros, J.A.O.; Sena-Cruz, J.M. A finite element model with discrete embedded elements for fibre reinforced composites. Comput. Struct. 2012,94–95, 22–33. [CrossRef] 10. Park, K.; Paulino, G.H.; Roesler, J. Cohesive fracture model for functionally graded fiber reinforced concrete. Cem. Conc. Res. 2010,40, 956–965. [CrossRef] 11. Yu, R.C.; Cifuentes, H.; Rivero, I.; Ruiz, G.; Zhang, X. Dynamic behaviour in fibre-reinforced cementitious composites. J. Mech. Phys. Solids 2016,93, 135–152. [CrossRef] 12. Radtke, F.K.F.; Simone, A.; Sluys, L.J. A partition of unity finite element method for simulating non-linear debonding and matrix failure in thin fibre composites. Int. J. Numer. Methods Eng. 2011 ,86, 453–476. [CrossRef] 13. Oliver, J.; Mora, D.F.; Huespe, A.E.; Weyler, R. A micromorphic model for steel fiber reinforced concrete. Int. J. Solids Struct. 2012,49, 2990–3007. [CrossRef] [PubMed] Materials 2017,10, 207 18 of 19 14. Montero-Chacón, F.; Schlangen, E.; Cifuentes, H.; Medina, F. A numerical approach for the design of multiscale fibre-reinforced cementitious composites. Phil. Mag. 2015,95, 3305–3327. [CrossRef] 15. Cundall, P.A. A computer model for simulating progressive large scale movements in blocky rock systems. Proc. Int. Symp. Rock Fract. 1971,1, 8–11. 16. Kawai, T. New discrete models and their application to seismic response analysis of structures. Nucl. Eng. Des. 1978,48, 207–229. [CrossRef] 17. Bolander, J.E.; Saito, S. Fracture analysis using spring network with random geometry. Eng. Fract. Mech. 1998,113, 1619–1630. 18. Zubelewicz, A.; Bažant, Z.P. Interface element modeling of fracture in aggregate composites. J. Eng. Mech.-ASCE 1987,113, 1619–1630. [CrossRef] 19. Bažant, Z.P.; Tabbara, M.R.; Kazemi, M.T.; Pijaudier-Cabot, G. Random particle model for fracture of aggregate or fiber composites. J. Eng. Mech.-ASCE 1990,116, 1686–1705. [CrossRef] 20. Cusatis, G.; Bažant, Z.P.; Cedolin, L. Confinement-shear lattice model for concrete damage in tension and compression: I. Theory. J. Eng. Mech.-ASCE 2003,129, 1439–1448. [CrossRef] 21. Schlangen, E.; Garboczi, E.J. Fracture simulations of concrete using lattice models: Computational aspects. Eng. Fract. Mech. 1997,57, 319–332. [CrossRef] 22. Bolander, J.E.; Saito, S. Discrete modeling of short-fibers reinforcement in cementitious composites. Adv. Cem. Mater. 1997,6, 76–86. [CrossRef] 23. Kunieda, M.; Ogura, H.; Ueda, N.; Nakamura, H. Tensile fracture process of Strain Hardening Cementitious Composites by means of three-dimensional meso-scale analysis. Cem. Concr. Comp. 2011 ,33, 956–965. [CrossRef] 24. Kang, J.; Kim, K.; Lim, Y.M.; Bolander, J.E. Modeling of fiber-reinforced cement composites: Discrete representation of fiber pullout. Int. J. Solids Struct. 2014,10, 1970–1979. [CrossRef] 25. Kozicki, J.; Tejchman, J. Effect of steel fibres on concrete behavior in 2D and 3D simulations using lattice model. Arch. Mech. 2010,62, 465–492. 26. Montero, F.; Schlangen, E. Modelling of fracture in fibre-cement based materials. Brittle Matrix Compos. 2012 , 10, 51–60. 27. Schauffert, E.; Cusatis, G. Lattice discrete particle model for fiber-reinforced concrete. I: Theory. J. Eng. Mech.-ASCE 2011,138, 826–833. [CrossRef] 28. Montero-Chacón, F.; Schlangen, E.; Medina, F. A lattice-particle approach for the simulation of fracture processes in fiber-reinforced high-performance concrete. In Proceedings of the VIII International Conference on Fracture Mechanics of Concrete and Concrete Structures, Toledo, Spain, 10–14 March 2013. 29. Montero-Chacón, F.; Medina, F. A lattice-particle approach to determine the RVE size for quasi-brittle materials. Eng. Comp. 2013,30, 246–262. [CrossRef] 30. Gitman, I.M.; Askes, H.; Sluys, L.J. Representative volume: Existence and size determination. Eng. Fract. Mech. 2007,74, 2518–2534. [CrossRef] 31. Nguyen, V.P.; Lloberas-Valls, O.; Stroeven, M.; Sluys, L.J. On the existence of representative volumes for softening quasi-brittle materials—A failure zone averaging scheme. Comp. Methods Appl. Mech. Eng. 2010 , 199, 3028–3038. [CrossRef] 32. Van Mier, J.G.M. Fracture Processes of Concrete, 1st ed.; CRC Press, Inc.: Boca Ratón, FL, USA, 1997. 33. Wriggers, P.; Moftah, S.O. Mesoscale models for concrete: Homogenisation and damage behavior. Finite Elem. Anal. Des. 2006,42, 623–636. [CrossRef] 34. Deeb, R.; Kulasegaram, S.; Karihaloo, B.L. 3D modelling of the flow of self-compacting concrete with or without steel fibres. Part I: Slump flow test. Comp. Part. Mech. 2014,1, 373–389. [CrossRef] 35. Ghebrab, T.T.; Soroushian, P. Mechanical properties of cement mortar: Development of structure-property relationships. Int. J. Concr. Struct. Mater. 2011,5, 3–10. [CrossRef] 36. Hsu, T.T.C.; Slate, F.O.; Sturman, G.M.; Winter, G. Microcracking of plain concrete and the shape of the stress-strain curve. J. ACI 1963,60, 209–224. 37. Rots, J.G.; Invernizzi, S. Regularized sequentially linear saw-tooth softening model. Int. J. Numer. Anal. Methods Geomech. 2004,28, 821–856. [CrossRef] 38. Kim, S.-M.; Abu Al-Rub, R.K. Meso-scale computational modeling of the plastic-damage response of cementitious composites. Cem. Concr. Res. 2011,41, 339–357. [CrossRef] Materials 2017,10, 207 19 of 19 39. Schellart, W.P. Shear test results for cohesion and friction coefficients for different granular materials: Scaling implications for their usage in analogue modeling. Tectonophysics 2000,324, 1–16. [CrossRef] 40. Nguyen, V.P.; Stroeven, M.; Sluys, L.J. Multiscale continuous and discontinuous modeling of heterogeneous materials: A review on recent developments. J. Multiscale Mod. 2011,3, 1–42. [CrossRef] 41. Cunha, V.; Barros, J.; Sena-Cruz, J. Pullout behavior of steel fibers in self-compacting concrete. J. Mater. Civ. Eng. 2010,22, 1–9. [CrossRef] 42. Kim, J.; Kim, D.; Kang, S.; Lee, J. Influence of sand to coarse aggregate ratio on the interfacial bond strength of steel fibers in concrete for nuclear power plant. Nucl. Eng. Des. 2012,252, 1–10. [CrossRef] 43. Karihaloo, B.L.; Lange-Kornbak, D. Optimization techniques for the design of high-performance fibre-reinforced concrete. Struct. Multidiscip. Optim. 2001,21, 32–39. [CrossRef] 44. Gopalaratnam, V.; Shah, S.P. Tensile failure of steel-fiber reinforced mortar. J. Eng. Mech.-ASCE 1987 ,113, 635–652. [CrossRef] 45. Johnston, C.D.; Coleman, R.A. Strength and deformation of steel fiber reinforced mortar in uniaxial tension. In Fiber Reinforced Concrete; ACI SP-44; American Concrete Institute: Farmington Hills, MI, USA, 1974; pp. 177–193. 46. ComitéEuro-International du Béton. CEB-FIP Model Code 1990, Bulletin D’Information; No. 213/214; Thomas Telford Services Ltd.: Lausanne, Switzerland, 1993. 47. Laranjeira, F.; Grünewald, S.; Walraven, J.; Blom, C.; Molins, C.; Aguado, A. Characterization of the orientation profile of steel fiber reinforced concrete. Mater. Struct. 2011,44, 1093–1111. [CrossRef] 48. SIMULIA Corp. Abaqus Theory Manual, version 6.8; Dassault Systémes: Providence, RI, USA, 2008. © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).