Numerical Simulation of the Fracture Behavior of High-Performance Fiber-Reinforced Concrete by Using a Cohesive Crack-Based Inverse Analysis

Fiber-reinforced concrete (FRC) has become an alternative for structural applications due its outstanding mechanical properties. The appearance of new types of fibres and the fibre cocktails that can be configured by mixing them has created FRC that clearly exceeds the minimum mechanical properties required in the standards. Consequently, in order to take full advantage of the contribution of the fibres in construction projects, it is of interest to have constitutive models that simulate the behaviour of the materials. This study aimed to simulate the fracture behaviour of five types of FRC, three with steel fibres, one with a combination of two types of steel fibers, and one with a combination of polyolefin fibres and two types of steel fibres, by means of an inverse analysis based on the cohesive crack approach. The results of the numerical simulations defined the softening functions of each FRC formulation and have pointed out the synergies that are created through use of fibre cocktails. The information supplied can be of help to engineers in designing structures with high-performance FRC.


Introduction
Fiber-reinforced concrete (FRC) can be considered a composite material as it merges two materials seeking to improve the properties that they individually have. While the plain concrete matrix provides the stiffness and the compressive strength, the fibres enhance the reduced tensile strength and brittleness of plain concrete. The result of this combination is a material with a remarkable compressive strength and an improved flexural and tensile strength, ductility, toughness, and impact resistance among other characteristics. [1][2][3]. Several studies have pointed out that the improvement of properties that the fibres provide to FRC depend on several factors such as their modulus of elasticity, constituent material, strength, fiber geometry, orientation, and dosage, among others [4][5][6].
Among all the available types of fibres suitable for being added to concrete, steel fibres, when added to concrete forming steel fiber-reinforced concrete (SFRC), have been the most used ones in structural applications such as footbridges, pavements, or tunnel linings [7][8][9][10][11]. In these applications, steel fibres prevent crack propagation and brittle failure, providing concrete with an improved strength and ductility after cracking. Moreover, polymer fibres such as PVA or polyolefin ones have proven to be a suitable alternative for structural uses and also might have fewer durability issues in chemically hazardous environments such as marine ones [12][13][14]. Another possibility that has been recently explored is using combinations of fibres, the so-called fibre cocktails, in order to obtain high performance FRC with tailor made properties. Several authors have studied the mechanical properties of fibre cocktails of micro and macro steel fibres with outstanding results [15][16][17]. Fibre cocktails composed by polymer and steel fibres have been also studied, being already available in papers dealing with concrete formulations that combined polypropylene and steel fibres. In most cases, such mixture aimed at complementing the structural characteristics of steel fibres with reduced drying and plastic shrinkage or even improved fire resistance [18,19].
As was mentioned in the previous paragraph, outstanding FRC properties might be obtained by mixing different types of fibers, but, in all cases, if a structural application is sought, the material properties must comply with certain requirements set in standards. There are national standards and structural design recommendations that establish the minimum mechanical properties that a certain type of FRC should meet to take advantage of the contribution of the fibres to the material properties [20][21][22]. However, in order to fully exploit the mechanical properties that the fiber cocktails provide, it is also necessary to supply reliable constitutive models that can be used in the structural design.
Due to the difficulties of developing such models, different theories have tried to describe the fracture process of concrete and FRC. Among them, the crack band theory [23] assumes that cracking takes place in a band where the cracks are homogeneously distributed, with the behaviour of the material described by a stress-strain relationship. Other authors have employed fracture mechanics to obtain the crack tip intensity factor by combining the fracture toughness of plain concrete and the stress-crack relationship provided by the fibre contribution [24]. Another possibility was to idealise the cracked material as a discontinuous entity where the crack that is described by a stress-crack opening relationship divides two sides of undamaged material that behaves elastically [25]. This approach has been employed with successful results in various published research. For instance, it was shown to be adequate not only for reproducing the fracture behavior of ultra-high-performance fibre-reinforced concrete notched prismatic specimens of several sizes reinforced with steel fibres [26] but also to determine the fracture energy of SFRC notched cylindrical specimens [27]. Moreover, the use of the cohesive crack approach has been employed for separating the contribution of the fracture energy and the fibre bridging zone in FRC [28]. In this study it was revealed that layered FRC might be more efficient for obtaining optimised structural elements. As a common feature of all these studies, it has been observed that the key factor of the model was the softening function that relates the stress-crack width relationship of the material.
To obtain such a relationship, some authors employed a numerical inverse approach that identified the mechanical and fracture properties required for performing an advanced analysis [29]. Other authors performed an analytical inverse analysis [30]. In all cases it was shown that the stress-crack width relationship that defines the behaviour of the cracks that appear in FRC depends on several parameters such as the plain concrete mechanical properties, fibre type [30], fibre dosage [31], and specimen size [32], among others.
Following that rationale, this paper sought to numerically replicate the three-point bending fracture behaviour of the same plain concrete with 0.89% volume fraction of OL and 3D and 5D steel fibers, respectively, by using multi-linear softening functions (OL fibres are short, straight, high-strength steel fibres, 3D fibres are steel fibers with hooked ends, and 5D fibers are high-strength steel fibres with double-hooked ends). To achieve this goal, experimental data found in literature were employed. The volume fraction was chosen, according to [33], in order to manufacture high-performance fiber-reinforced concrete mixes able to broadly exceed the requirements set in [20][21][22]. To simulate the material behaviour, the suitability of using multilinear softening functions to represent the σ-w relation of the mentioned types of SFRC had to be determined. Moreover, once the previous tasks were fulfilled, the constitutive model of the same plain concrete with a binary and ternary cocktail of fibres was found. In addition, a discussion regarding the differences among the softening functions and their relationship with the fibre geometry was carried out.
The significance of this research relies on the potential optimization of the structural design that might be achieved by implementing the constitutive models in structural design software. Moreover, once the individual softening functions of each SFRC have been determined, the analytical deduction of the softening functions of the cocktails of fibres may be explored in future contributions. Another point that should be underlined is that the mentioned optimisation in the use of the material properties might result in remarkable reductions, or even a total substitution, of the steel-bar reinforcement without compromising the reliability of FRC structural elements. Besides, the reduction of the amount of steel used could suppose not only a decrease in the economic costs but also in an increment of the sustainability of the structure due to the savings in raw materials and safety and improvements of durability.

Material Production and Experimental Campaign
This section will briefly present the manufacturing process of the concrete formulations, the fracture behavior of the mixes was simulated in the following sections. A detailed description of the materials used and the steps taken can be seen in [33].
Five formulations of self-compacting FRC were prepared, adding several fiber combinations to the same plain concrete formulation. Such plain concrete was designed considering three aggregate types: two types of gravel with particle sizes between 4-8 mm and 4-12 mm, respectively, and one type of sand with particle sizes between 0-2 mm. The maximum size of the aggregates was 12.7 mm. In addition, limestone powder was employed, seeking to enhance the workability of the mixes. In all formulations a Portland cement-type EN 197-1 CEM I 52.5 was used as binder. A policarboxylic superplasticiser was added to the formulations pursuing to achieve the self-compactability of all types of concrete while in a fresh state. Tap water was also employed in concrete production. The proportions of all components can be seen in Table 1. Three types of fibers (OL, 3D and 5D) were used in order to obtain three types of SFRC. Such formulations were termed SFRC-OL, SFRC-3D and SFRC-5D, respectively. In addition, two fibre cocktails were employed in the other two formulations. One of the hybrid formulations, which was named H1, was obtained by mixing OL and 5D fibres. The other hybrid formulation employed OL, 5D, and polyolefin fibers (PF). This last formulation was called H2. The fibers' characteristics and their appearance can be seen in Table 2 and Figure 1, respectively.  The five concrete formulations were designed seeking to obtain high-performance concrete that was capable of surpassing the requirements set in the most common standards [20][21][22] for reducing or even eliminating the steel reinforcing bars. Consequently, the volume fraction of fibres added was equal to 0.89% in all concrete types. It should be pointed out that all the formulations were designed to meet the requirements that are set for self-compacting concrete while in a fresh state. As it was previously stated, three concrete types were SFRC with only one type of fibre and the other two were concrete types where a cocktail of fibres had been added. In H1 formulation, the volume fraction of OL and 5D fibres were equal to 0.445% and their sum was equal to 0.89%. Similarly, H2 had equal volume fractions of OL, 5D and PF, with each one being 0.293%, which is a volume fraction of 0.89% combined. It was observed that the experimental campaign had as a main target to characterize the mechanical behaviour of high-performance FRC and also analyse the synergies and improvements that could be obtained by combining several fibre types.
All concrete batches were manufactured in a planetary mixer with vertical axis following the same scheme. A thorough description of the process can be seen in [33]. From each batch, six cylindrical specimens of 150 mm of diameter and 300 mm of height were cast. In addition, four prismatic specimens of dimensions and 600 × 150 × 150 mm³ were also manufactured. All specimens were cured in a climatic chamber until 28 days of age. The cylindrical specimens were used to obtain the compressive strength, the modulus of elasticity, and the tensile strength following EN 12390-3 [34], EN 12390-13 [35], and EN 12390-6 [36], respectively. All tests were conducted at 28 days of age of the material. In Table 3, the results of the mentioned tests can be seen. As can be seen in Table 3, there were no remarkable differences in the compressive strength of all the formulations tested. Likewise, no notable changes of the tensile strength were found among the mixes. On the contrary, regarding the modulus of elasticity, it should be highlighted that H1 boasted a modulus of elasticity that was clearly superior to The five concrete formulations were designed seeking to obtain high-performance concrete that was capable of surpassing the requirements set in the most common standards [20][21][22] for reducing or even eliminating the steel reinforcing bars. Consequently, the volume fraction of fibres added was equal to 0.89% in all concrete types. It should be pointed out that all the formulations were designed to meet the requirements that are set for self-compacting concrete while in a fresh state. As it was previously stated, three concrete types were SFRC with only one type of fibre and the other two were concrete types where a cocktail of fibres had been added. In H1 formulation, the volume fraction of OL and 5D fibres were equal to 0.445% and their sum was equal to 0.89%. Similarly, H2 had equal volume fractions of OL, 5D and PF, with each one being 0.293%, which is a volume fraction of 0.89% combined. It was observed that the experimental campaign had as a main target to characterize the mechanical behaviour of high-performance FRC and also analyse the synergies and improvements that could be obtained by combining several fibre types.
All concrete batches were manufactured in a planetary mixer with vertical axis following the same scheme. A thorough description of the process can be seen in [33]. From each batch, six cylindrical specimens of 150 mm of diameter and 300 mm of height were cast. In addition, four prismatic specimens of dimensions and 600 × 150 × 150 mm 3 were also manufactured. All specimens were cured in a climatic chamber until 28 days of age. The cylindrical specimens were used to obtain the compressive strength, the modulus of elasticity, and the tensile strength following EN 12390-3 [34], EN 12390-13 [35], and EN 12390-6 [36], respectively. All tests were conducted at 28 days of age of the material. In Table 3, the results of the mentioned tests can be seen. As can be seen in Table 3, there were no remarkable differences in the compressive strength of all the formulations tested. Likewise, no notable changes of the tensile strength were found among the mixes. On the contrary, regarding the modulus of elasticity, it should be highlighted that H1 boasted a modulus of elasticity that was clearly superior to the ones found in the other concrete types. However, for the purpose of this manuscript, such changes were of minor influence, as the modulus of elasticity only possessed a remarkable influence in the initial part of the fracture behaviour under a three-point bending setup.
The fracture behaviour was determined following the process shown in [37]. The results of the tests boasted a limited scattering, and the average of the results can be seen in Figure 2.
such changes were of minor influence, as the modulus of elasticity only possessed a remarkable influence in the initial part of the fracture behaviour under a three-point bending setup.
The fracture behaviour was determined following the process shown in [37]. The results of the tests boasted a limited scattering, and the average of the results can be seen in Figure 2. In Figure 2a it can be seen that, when the fiber type changed, the fracture behavior of the SFRC studied significantly varied. It can be seen that the overall best performance was obtained with the 3D fibres, followed by the 5D fibres and lastly by the OL fibres. SFRC-3D reached the highest load (maximum load and remaining load) followed by SFRC-5D and SFRC-OL. It is worth noting that, in SFRC-5D and SFRC-3D, alike, the maximum load appeared around 1 mm of CMOD, while in SFRC-OL, the maximum load was reached at a CMOD around 0.125 mm. After reaching the maximum load, all concrete specimens showed a progressive unloading.
The fracture behaviour of the concrete mixes where cocktail of fibres were added can be seen in Figure 2b.
In Figure 2a,b the curves shown correspond to at least the average behaviour of three specimens. It should be also highlighted that no remarkable scattering was found in the tests. Regarding the fracture behaviour of the mentioned formulations, it can be stated that the cocktails used in H1 and H2 clearly enhanced the fracture behaviour of any of the formulations that had only one type of fibres.
In order to determine if it is possible to consider the contribution of the fibres in the structural design and, consequently reduce the amount of steel reinforcement, [22] sets of the strength of the material at a CMOD value of 0.5 mm, which is often known as fR1, should be at least 40% of the load at the limit of proportionality (fLOP). Moreover, it should be also checked that the load at a CMOD value of 2.5 mm should be at least 20% of fLOP, usually termed fR3. In Table 4, it can be observed that all concrete formulations exceeded the requirements set by [22].  In Figure 2a it can be seen that, when the fiber type changed, the fracture behavior of the SFRC studied significantly varied. It can be seen that the overall best performance was obtained with the 3D fibres, followed by the 5D fibres and lastly by the OL fibres. SFRC-3D reached the highest load (maximum load and remaining load) followed by SFRC-5D and SFRC-OL. It is worth noting that, in SFRC-5D and SFRC-3D, alike, the maximum load appeared around 1 mm of CMOD, while in SFRC-OL, the maximum load was reached at a CMOD around 0.125 mm. After reaching the maximum load, all concrete specimens showed a progressive unloading.
The fracture behaviour of the concrete mixes where cocktail of fibres were added can be seen in Figure 2b.
In Figure 2a,b the curves shown correspond to at least the average behaviour of three specimens. It should be also highlighted that no remarkable scattering was found in the tests. Regarding the fracture behaviour of the mentioned formulations, it can be stated that the cocktails used in H1 and H2 clearly enhanced the fracture behaviour of any of the formulations that had only one type of fibres.
In order to determine if it is possible to consider the contribution of the fibres in the structural design and, consequently reduce the amount of steel reinforcement, [22] sets of the strength of the material at a CMOD value of 0.5 mm, which is often known as f R1, should be at least 40% of the load at the limit of proportionality (f LOP ). Moreover, it should be also checked that the load at a CMOD value of 2.5 mm should be at least 20% of f LOP , usually termed f R3 . In Table 4, it can be observed that all concrete formulations exceeded the requirements set by [22]. Table 4. Residual strength of the concrete formulation in relation with the requirements set by [22]. Similarly to what was mentioned in the case of the production process, a detailed description and discussion of the results can be found in [33].

Numerical Simulations
The concept of the cohesive crack zone [25,38,39] was implemented in an element formulation, with an embedded cohesive crack, in order to reproduce the cracking process of the concrete formulations. Although this concept was firstly developed for plain concrete [40], it has been subsequently applied in fibre cementitious materials with successful results [41].
This approach considers that fracture takes place in mode I, and, hence, it was assumed that the cohesive stress vector t transmitted across the crack faces was parallel to the crack displacement vector w (central forces model). If the magnitude of the crack opening vector |w| does not decrease, the relation between both parameters can be stated as appears in (1) where f (|w|) is the softening function for pure opening, which, in this manuscript, will be defined by multilinear functions, where unloading and reloading branches are aligned with the origin and the softening function is defined by a set of points.
The mentioned formulation is referred to as constant strain triangular finite elements that only have one integration point. It should be mentioned that only the directions parallel to the triangle sides are considered as cracking directions.
If the element is cracked and, consequently, the crack direction is determined, the stress vector t can be obtained as It should be mentioned that the stress vector t is constant along the crack, with h being the triangle height, A representing the area of the element, L representing the crack length, σ representing the stress tensor, and n representing the unit vector normal to that side and to the crack. A more detailed explanation of the implementation can be found in [41].
In the cohesive crack theory, it is assumed that outside the crack the material remains elastic. The stress tensor can be obtained using Equation (3), where the inelastic behaviour is subtracted to correct the elastic deformation of the element by including the crack displacement.
In Equation (3), E stands for the elastic tangent tensor, ε α is the apparent strain vector obtained with the nodal displacements, and b + is the gradient vector corresponding to the shape function of the solitary node, as appears in Equation (4): If it is considered that t = σ n and also using Equations (1) and (3), the following expression could be obtained.
where 1 stands for the identity tensor. If an iterative algorithm is used, w could be found to satisfy Equation (5). The approach previously explained was introduced in a material subroutine in a commercial finite element code (ABAQUS) seeking to reproduce the material behaviour. An external file with the coordinates of the nodes and the elements is required to feed the subroutine with the geometry of the elements. If the material has not reached the tensile strength (f ct ), the element behaves elastically. The tensile strength influences the limit of proportionality of the experimental curves, being greater as f ct increases. Nevertheless, if such value is exceeded, the direction of the maximum principal stress is found and a crack perpendicular to such direction is introduced. From that point onwards, the behaviour of the element is governed, if strains keep growing, by the softening function of the material, which relates the stress and the crack opening relation (σ-w). Apart from the tensile strength, the subroutine requires the modulus of elasticity of the material E. Although it is true that the presence of fibres in concrete changes the post-peak behaviour of the material when subjected to compressive stresses, given that for high crack widths the material behaviour is mainly under tensile stresses, no compressive damage was considered in the material implementation. Moreover, no experimental data appeared in [33] regarding the post-peak compressive strength of any of the mixes analysed.
The shape and characteristic points of the softening functions of all the mixes analysed were determined by the authors in the current study based on an iterative inverse analysis. The best reproduction of the fracture test results was sought. Such process can be summarised as the following steps. Firstly, the softening function is implemented in the element and the points that define it are set. Secondly, the simulation is performed and the behaviour of the numerical model when subjected to the fracture test is obtained. Thirdly, the experimental and the numerical behaviors are compared. This can be carried out using the load-deflection or the load-crack mouth opening displacement (CMOD) curves. If the numerical behavior does not reproduce the experimental behaviour, the points that define the softening function are modified and the next iteration starts again. If the numerical behaviour simulates with a certain degree of accuracy the experimental behavior, the points proposed in the first step are considered valid.
In order to reproduce the material behavior, a plane mesh with plane stress conditions was employed. The mesh was composed of constant strain triangles, which are finer in the whereabouts of the notch. An image of one of the meshes can be seen in Figure 3. This scheme of mesh was considered valid in [32,42].
subroutine with the geometry of the elements. If the material has not reached the tensile strength (fct), the element behaves elastically. The tensile strength influences the limit of proportionality of the experimental curves, being greater as fct increases. Nevertheless, if such value is exceeded, the direction of the maximum principal stress is found and a crack perpendicular to such direction is introduced. From that point onwards, the behaviour of the element is governed, if strains keep growing, by the softening function of the material, which relates the stress and the crack opening relation (σ-w). Apart from the tensile strength, the subroutine requires the modulus of elasticity of the material E. Although it is true that the presence of fibres in concrete changes the post-peak behaviour of the material when subjected to compressive stresses, given that for high crack widths the material behaviour is mainly under tensile stresses, no compressive damage was considered in the material implementation. Moreover, no experimental data appeared in [33] regarding the post-peak compressive strength of any of the mixes analysed.
The shape and characteristic points of the softening functions of all the mixes analysed were determined by the authors in the current study based on an iterative inverse analysis. The best reproduction of the fracture test results was sought. Such process can be summarised as the following steps. Firstly, the softening function is implemented in the element and the points that define it are set. Secondly, the simulation is performed and the behaviour of the numerical model when subjected to the fracture test is obtained. Thirdly, the experimental and the numerical behaviors are compared. This can be carried out using the load-deflection or the load-crack mouth opening displacement (CMOD) curves. If the numerical behavior does not reproduce the experimental behaviour, the points that define the softening function are modified and the next iteration starts again. If the numerical behaviour simulates with a certain degree of accuracy the experimental behavior, the points proposed in the first step are considered valid.
In order to reproduce the material behavior, a plane mesh with plane stress conditions was employed. The mesh was composed of constant strain triangles, which are finer in the whereabouts of the notch. An image of one of the meshes can be seen in Figure 3. This scheme of mesh was considered valid in [32,42].

Results and Discussion
In Figure 4a, the comparison of the results obtained in the experimental campaign and the simulations can be seen regarding SFRC-OL. As it can be clearly observed, the degree of accuracy that was achieved was remarkable, with the only slight difference present in the peak load. In contrast to what was shown in [41], the softening function in this occasion was composed by four straight stretches. The first one started when the opening of the crack was null and the stress was equal to the tensile strength registered in tests

Results and Discussion
In Figure 4a, the comparison of the results obtained in the experimental campaign and the simulations can be seen regarding SFRC-OL. As it can be clearly observed, the degree of accuracy that was achieved was remarkable, with the only slight difference present in the peak load. In contrast to what was shown in [41], the softening function in this occasion was composed by four straight stretches. The first one started when the opening of the crack was null and the stress was equal to the tensile strength registered in tests carried out in plain concrete. There was a sudden drop of the load-bearing capacity of the material in the initial stretch that might have been related to the crack width required by the fibres to contribute to the strength of SFRC-OL. Once the fibres started to enhance the mechanical properties of the concrete matrix, the load-bearing capacity decrease was more gradual. Such a feature was also perceived in the last two stretches of the softening function.
carried out in plain concrete. There was a sudden drop of the load-bearing capacity of the material in the initial stretch that might have been related to the crack width required by the fibres to contribute to the strength of SFRC-OL. Once the fibres started to enhance the mechanical properties of the concrete matrix, the load-bearing capacity decrease was more gradual. Such a feature was also perceived in the last two stretches of the softening function. Similarly to Figure 4, in Figure 5a we can see the comparison between the experimental results registered when testing SFRC-3D in a three-point bending fracture test and the results that were obtained in the numerical simulations. It should be highlighted that the accuracy of the simulation was remarkable not only in the peak load but also throughout the unloading branch of the test. Moreover, even at deflection values beyond 6 mm, the simulation was able to replicate the behaviour of SFRC-3D with precision. Regarding the softening function that enabled such accuracy, which is shown in Figure 5b, it should be underlined that, in contrast to what was shown in Figure 4b, there was no initial steep unloading for reduced crack widths. On the contrary, as can be perceived in Figure 5b, the initial part of the softening function maintained the load-bearing capacity close to the tensile strength of plain concrete. Such a feature might be related with the shape of the fibres that boasted hooked ends. Once this initial stretch was surpassed, there was a progressive loss in the load-bearing capacity of SFRC-3D, which is notoriously steeper than the first stretch. Similarly to Figure 4, in Figure 5a we can see the comparison between the experimental results registered when testing SFRC-3D in a three-point bending fracture test and the results that were obtained in the numerical simulations. It should be highlighted that the accuracy of the simulation was remarkable not only in the peak load but also throughout the unloading branch of the test. Moreover, even at deflection values beyond 6 mm, the simulation was able to replicate the behaviour of SFRC-3D with precision. Regarding the softening function that enabled such accuracy, which is shown in Figure 5b, it should be underlined that, in contrast to what was shown in Figure 4b, there was no initial steep unloading for reduced crack widths. On the contrary, as can be perceived in Figure 5b, the initial part of the softening function maintained the load-bearing capacity close to the tensile strength of plain concrete. Such a feature might be related with the shape of the fibres that boasted hooked ends. Once this initial stretch was surpassed, there was a progressive loss in the load-bearing capacity of SFRC-3D, which is notoriously steeper than the first stretch.
carried out in plain concrete. There was a sudden drop of the load-bearing capacity of the material in the initial stretch that might have been related to the crack width required by the fibres to contribute to the strength of SFRC-OL. Once the fibres started to enhance the mechanical properties of the concrete matrix, the load-bearing capacity decrease was more gradual. Such a feature was also perceived in the last two stretches of the softening function. Similarly to Figure 4, in Figure 5a we can see the comparison between the experimental results registered when testing SFRC-3D in a three-point bending fracture test and the results that were obtained in the numerical simulations. It should be highlighted that the accuracy of the simulation was remarkable not only in the peak load but also throughout the unloading branch of the test. Moreover, even at deflection values beyond 6 mm, the simulation was able to replicate the behaviour of SFRC-3D with precision. Regarding the softening function that enabled such accuracy, which is shown in Figure 5b, it should be underlined that, in contrast to what was shown in Figure 4b, there was no initial steep unloading for reduced crack widths. On the contrary, as can be perceived in Figure 5b, the initial part of the softening function maintained the load-bearing capacity close to the tensile strength of plain concrete. Such a feature might be related with the shape of the fibres that boasted hooked ends. Once this initial stretch was surpassed, there was a progressive loss in the load-bearing capacity of SFRC-3D, which is notoriously steeper than the first stretch. In Figure 6a, we can observe the comparison between the experimental results of the three-point bending fracture tests of SFRC-5D and their numerical reproduction. The dashed line that represents the numerical simulation clearly follows the experimental trends with accuracy. The peak load was captured with precision and even the progressive loss of stiffness that the material suffered before reaching such point was reproduced. Moreover, there was an almost identical unloading in the simulation and in the experimental test from the peak load to the end of the test. Concerning the softening function that enabled such an accurate simulation, it should be mentioned that, in contrast with what was shown in Figure 5b, the softening function of SFRC-5D that can be seen in Figure 6b boasts a sudden decrease of the load-bearing capacity after surpassing the tensile strength. This phenomenon was not expected because such a load drop did not appear in the case of 3D fibres. If the crack width continued to grow, the load-bearing capacity would increase up to a local maximum, which would be followed by the final unloading of the material.
three-point bending fracture tests of SFRC-5D and their numerical reproduction. The dashed line that represents the numerical simulation clearly follows the experimental trends with accuracy. The peak load was captured with precision and even the progressive loss of stiffness that the material suffered before reaching such point was reproduced. Moreover, there was an almost identical unloading in the simulation and in the experimental test from the peak load to the end of the test. Concerning the softening function that enabled such an accurate simulation, it should be mentioned that, in contrast with what was shown in Figure 5b, the softening function of SFRC-5D that can be seen in Figure  6b boasts a sudden decrease of the load-bearing capacity after surpassing the tensile strength. This phenomenon was not expected because such a load drop did not appear in the case of 3D fibres. If the crack width continued to grow, the load-bearing capacity would increase up to a local maximum, which would be followed by the final unloading of the material. In this work, not only the experimental curves of the three-point bending fracture tests of the formulations with the addition of one fiber were simulated but also the formulations with an addition of cocktail of fibres. In the case of the formulation H1, the cocktail of fibres was composed by OL fibres and 5D fibres. At this point, it should be mentioned that the fracture curve was simulated following the same iterative process that was applied to the concrete formulation with only one type of fibre. Consequently, although the softening functions of SFRC-OL and SFRC-5D were known, the softening function of H1 could not be inferred from them.
In Figure 7a, we can see the comparison between the three-point bending fracture tests of H1 and the results obtained in the simulations. It can be observed that the simulations with the softening function, shown in Figure 7b, posed a high resemblance with the experimental results, both in the loading and unloading part of the curve. Regarding the softening function, it should be underlined that it featured the initial sudden load-bearing capacity drop typical of 5D fibres. Besides, the same steep recovery appeared right after the initial unloading took place. After reaching the maximum of this second stretch, there was a progressive unloading until the total collapse of the material. It should be also mentioned that the maximum crack width in the softening function of H1 was larger than any of the concrete formulations that had only one type of fibre. In this work, not only the experimental curves of the three-point bending fracture tests of the formulations with the addition of one fiber were simulated but also the formulations with an addition of cocktail of fibres. In the case of the formulation H1, the cocktail of fibres was composed by OL fibres and 5D fibres. At this point, it should be mentioned that the fracture curve was simulated following the same iterative process that was applied to the concrete formulation with only one type of fibre. Consequently, although the softening functions of SFRC-OL and SFRC-5D were known, the softening function of H1 could not be inferred from them.
In Figure 7a, we can see the comparison between the three-point bending fracture tests of H1 and the results obtained in the simulations. It can be observed that the simulations with the softening function, shown in Figure 7b, posed a high resemblance with the experimental results, both in the loading and unloading part of the curve. Regarding the softening function, it should be underlined that it featured the initial sudden load-bearing capacity drop typical of 5D fibres. Besides, the same steep recovery appeared right after the initial unloading took place. After reaching the maximum of this second stretch, there was a progressive unloading until the total collapse of the material. It should be also mentioned that the maximum crack width in the softening function of H1 was larger than any of the concrete formulations that had only one type of fibre. Regarding the fracture behaviour of H2, in Figure 8a we can observe the precision obtained in the simulations if such results are compared with the experimental ones. It should be emphasised that the numerical simulations were able to reproduce not only the initial loading stretch but also the progressive loss of stiffness that the material suffered after reaching 20 kN until the peak load of the test was registered. It should be also men- Regarding the fracture behaviour of H2, in Figure 8a we can observe the precision obtained in the simulations if such results are compared with the experimental ones. It should be emphasised that the numerical simulations were able to reproduce not only the initial loading stretch but also the progressive loss of stiffness that the material suffered after reaching 20 kN until the peak load of the test was registered. It should be also mentioned that, as happened in the SFRC-5D and SFRC-3D tests, the peak load was registered gradually in contrast to the abrupt loss of stiffness observed in SFRC-OL. Such a feature entails the presence of 5D fibers capable of sustaining loads for reduced deflections. These two characteristics of the experimental curve were notably well captured by the numerical simulations. Moreover, a great similarity of the experimental and the numerical curves were seen during the complete test. If Figure 8b is analyzed, it can be mentioned that there was a sharp decrease of the load-bearing capacity of the material right after surpassing the tensile strength of plain concrete. Nevertheless, it can be also said that the presence of 5D and PF fibres enabled H2 to partially recover the load-bearing capacity even at low values of crack openings. Considering the notorious differences between the modulus of elasticity of 5D and PF fibres, 5D fibres might be considered the main responsible for such recovery. Another characteristic that should be underlined is that the maximum crack opening for this cocktail of fibres was remarkably greater than the ones shown in SFRC-3D, SFRC-5D, and SFRC-OL formulations. Regarding the fracture behaviour of H2, in Figure 8a we can observe the precision obtained in the simulations if such results are compared with the experimental ones. It should be emphasised that the numerical simulations were able to reproduce not only the initial loading stretch but also the progressive loss of stiffness that the material suffered after reaching 20 kN until the peak load of the test was registered. It should be also mentioned that, as happened in the SFRC-5D and SFRC-3D tests, the peak load was registered gradually in contrast to the abrupt loss of stiffness observed in SFRC-OL. Such a feature entails the presence of 5D fibers capable of sustaining loads for reduced deflections. These two characteristics of the experimental curve were notably well captured by the numerical simulations. Moreover, a great similarity of the experimental and the numerical curves were seen during the complete test. If Figure 8b is analyzed, it can be mentioned that there was a sharp decrease of the load-bearing capacity of the material right after surpassing the tensile strength of plain concrete. Nevertheless, it can be also said that the presence of 5D and PF fibres enabled H2 to partially recover the load-bearing capacity even at low values of crack openings. Considering the notorious differences between the modulus of elasticity of 5D and PF fibres, 5D fibres might be considered the main responsible for such recovery. Another characteristic that should be underlined is that the maximum crack opening for this cocktail of fibres was remarkably greater than the ones shown in SFRC-3D, SFRC-5D, and SFRC-OL formulations.  The accuracy of the simulations performed in all the concrete mixes was assessed by comparing the fracture energies that were obtained in the experimental tests and the numerical simulations. In the last column of Table 5, the relative error of the fracture energy simulated can be seen, being in all cases below 2%.

Discussion
In Figure 9a, a comparison of the softening functions of the SFRC formulations with only one type of fibre addition can be observed. It is obvious that there were significant differences among the three curves. The first feature that can be noted is that, in the case of the SFRC-OL and SFRC-5D formulations, right after surpassing the tensile strength, a sudden drop of the load bearing capacity of the material appeared. While this might be expected in the case of the SFRC-OL formulation due to the straight shape of the OL fibres, in the case of the 5D fibers, the double hook present in the ends of the fibres might have prevented any unloading from happening. Moreover, when seeking to determine the softening function of SFRC-5D, it became evident that the slight reduction of the stiffness that appeared in the load-deflection curves shown in Figure 6a was greatly influenced by this stretch of the softening function. In contrast to what was previously explained, the first stretch of the SFRC-3D softening function showed that the fibres could bear the load that the concrete matrix transferred to them when the tensile strength of the concrete was reached. Thus, such a stretch only showed a slight unloading. In the case of the softening functions of H1 and H2 formulations, it can be observed in Figure 9b that they boasted similar shapes. Both of them had a rapid reduction of the load-bearing capacity right after reaching the tensile strength of the concrete matrix. This might be explained by the presence of 5D and OL fibres in both formulations. Moreover, in the case of H2, the presence of PF fibres in the fibre cocktail may have been also related with this feature. Regarding the re-loading stage of both softening functions, the H1 formulation was capable of gaining stress up to values close to the tensile strength of plain concrete. This feature might have been caused by the greater proportion of 5D fibres present in the H1 mix. On the contrary, H2 boasted a re-loading stage more subtle and with less slope than the corresponding stretch of H1. As was stated in [43], the softening functions that can represent the fracture behavior of polyolefin fibre-reinforced concrete (PFRC) mixes boasted a reloading stage after the first unloading. Therefore, such changes might have accounted for the presence of a greater number of 5D fibres in H1 and also for the lower modulus of elasticity of the PF fibres present in H2. Another feature that should be underlined is that the softening functions of H1 and H2 alike had maximum crack widths that were greater than any of the softening functions of the concretes, where only one type of fiber was added. This increment of the maximum crack width might be an effect of the synergies that appear when cocktails of fibres are used.
A comparison among the points that define the softening functions of the five mixes analysed can be seen in Table 6. The mentioned features can be clearly perceived comparing the coordinates of the points that define the softening functions. After the first stretch, the softening function of the SFRC-OL fibers continued to develop a progressive unloading until the failure of the material took place. It should be mentioned that the slope of the rest of the unloading stages was much more progressive than the first stretch. This feature conferred SFRC-OL a certain ductility. Regarding the remaining part of the SFRC-3D softening function, it can be seen in Figure 9a that there was a load drop stretch until the collapse of the material happens. On the contrary, the softening function of SFRC-5D boasted an increment of the load-bearing capacity in the second stretch of the curve. Such an increment might be a result of the anchorage between the fibres and the concrete matrix. Once the anchorage fails, the material begins to unload until the failure is produced. Another feature that should be mentioned is that, in all concrete formulations, the maximum crack width was smaller than the length of the fibres added, with no evident relation to the fiber length. Such a relation was evident in the case of OL fibers, the length of which was 13 mm and the maximum crack width if which was 5 mm. Similarly, the length of 5D fibres was equal to 60 mm, while the maximum crack width was 8.3 mm. An analogue comment can be performed regarding 3D fibres.
In the case of the softening functions of H1 and H2 formulations, it can be observed in Figure 9b that they boasted similar shapes. Both of them had a rapid reduction of the load-bearing capacity right after reaching the tensile strength of the concrete matrix. This might be explained by the presence of 5D and OL fibres in both formulations. Moreover, in the case of H2, the presence of PF fibres in the fibre cocktail may have been also related with this feature. Regarding the re-loading stage of both softening functions, the H1 formulation was capable of gaining stress up to values close to the tensile strength of plain concrete. This feature might have been caused by the greater proportion of 5D fibres present in the H1 mix. On the contrary, H2 boasted a re-loading stage more subtle and with less slope than the corresponding stretch of H1. As was stated in [43], the softening functions that can represent the fracture behavior of polyolefin fibre-reinforced concrete (PFRC) mixes boasted a reloading stage after the first unloading. Therefore, such changes might have accounted for the presence of a greater number of 5D fibres in H1 and also for the lower modulus of elasticity of the PF fibres present in H2. Another feature that should be underlined is that the softening functions of H1 and H2 alike had maximum crack widths that were greater than any of the softening functions of the concretes, where only one type of fiber was added. This increment of the maximum crack width might be an effect of the synergies that appear when cocktails of fibres are used.
A comparison among the points that define the softening functions of the five mixes analysed can be seen in Table 6. The mentioned features can be clearly perceived comparing the coordinates of the points that define the softening functions.

Conclusions
This study implemented multi-linear softening functions in a user material subroutine of a finite element code (ABAQUS) in order to simulate the three-point bending fracture behaviour of three SFRC and two FRC with cocktails of fibres. The results obtained showed that the methodology used is a robust and efficient approach to simulate the fracture behaviour of the concrete mixes analysed.
The mentioned approach reproduced the fracture behaviour of the five concrete types studied with accuracy, obtaining less than a 2% difference between the experimental and the simulated fracture energy.
The analysis of the softening function corresponding to SFRC-5D showed that the presence of the double-hooked end in the 5D fibres might not provide an increment of the fiber-matrix anchorage for reduced CMOD values. However, in the case of 3D fibers, no load drops appeared at initial crack widths of the softening function. Regarding SFRC-OL, its softening function boasted a progressive unloading of the material. Lastly, no relationship was observed between the length of the fibres and the maximum crack width of the softening functions found.
Regarding the softening functions of H1 and H2, the presence of OL and 5D fibres in both cocktails explains the presence of an unloading stretch and a re-loading stage. The synergies that appeared when combining several types of fibres generated an increment of the maximum crack width both in H1 and H2. However, a further analysis should be performed for a better measurement of the enhancement that the combination of fibers can provide.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.