Predicting the Tensile Behaviour of Ultra-High Performance Fibre-Reinforced Concrete from Single-Fibre Pull-Out Tests

In this paper, a prediction model for the tensile behaviour of ultra-high performance fibre-reinforced concrete is proposed. It is based on integrating force contributions of all fibres crossing the crack plane. Piecewise linear models for the force contributions depending on fibre orientation and embedded length are fitted to force–slip curves obtained in single-fibre pull-out tests. Fibre characteristics in the crack are analysed in a micro-computed tomography image of a concrete sample. For more general predictions, a stochastic fibre model with a one-parametric orientation distribution is introduced. Simple estimators for the orientation parameter are presented, which only require fibre orientations in the crack plane. Our prediction method is calibrated to fit experimental tensile curves.

The post-crack tensile strength of UHPFRC highly depends on geometric characteristics of the fibre system such as fibre content, shape, aspect ratio, spatial arrangement, and orientation. The latter two may vary depending on production parameters of the concrete [13]. A high tensile strength is obtained if the crack-crossing fibres are evenly distributed over the crack and if the fibres are aligned to the tensile axis [13][14][15].
Several authors [16][17][18][19][20][21][22][23][24][25][26][27] proposed prediction models for the tensile behaviour based on single-fibre pull-out force contributions. The prediction models by Li et al. [16] and Wuest et al. [17] integrate the force contributions of crack-crossing fibres over the entire crack plane in the tensioned fibre reinforced concrete. Such a sectional analysis provides a relation between stress and crack opening. The force contribution of individual crackcrossing fibres can be approximated from single-fibre pull-out tests for selected fibre orientations and embedded lengths [18,21,[27][28][29] or by analytical equations [19,20,23,24]. For modelling purposes, the experimental single-fibre pull-out curves are averaged [25] or idealised [23,27] as piecewise linear with cut points equal to the tensile force at the end of the linear phase and the ultimate force. Besides the fibres, also the matrix contributes to the tensile resistance. Experimental fibre-pull-out curves naturally contain information on fibrematrix bond properties. If such data are not available, the matrix resistance contribution can be modelled by bilinear [30,31], trilinear [32,33], or exponential [22] functions and added to the fibre force contributions.
A common assumption in calculating the force contributions of crack-crossing fibres is that the fibre orientation is uniformly distributed [19,22,31,34,35], i.e., that concrete is an isotropic material-this simplifies calculations. In practice, however, the fibre orientation distribution deviates significantly from this assumption [13]. This should be taken into account for formulating more realistic models. Empirical fibre orientation distributions determined by image analysis are used in [25,27,36] while [37] derive orientation information from hydrodynamic equations related to concrete flow. The authors are not aware of any study that models the fibre direction distribution stochastically beyond the previously mentioned simple uniformity assumptions or empirical investigations.
An alternative prediction approach is to derive simple analytical equations based on concrete characteristics. A review of such models for predicting the shear capacity is given by [38]. None of these take the fibre orientation distribution into account.
Furthermore, (nonlinear) finite element approaches are nowadays popular, for instance, to predict the overall hysteretic response of cyclic-loaded SFRC [39] or to verify the micromechanics-based fibre-bridging curves of tensioned UHPFRC [23]. Drawbacks of FE methods are the higher computational effort compared to analytical methods and that values of material parameters (such as friction parameters) have to be available. Simulating the contact relationship between steel fibre and matrix via contact [40,41] and coupling [42,43] algorithms requires identical locations for fibre and matrix nodes. This is impractical in the mesh generation and for performing numerical calculations incorporating large quantities of small fibres [35]. The recently introduced model of [35] overcomes this problem, but again only considers uniformly distributed fibre orientations.
Combining suitable imaging techniques with quantitative image analysis is an established way to measure geometric characteristics of the fibre system. For the prediction of tensile behaviour, characteristics of fibres intersecting the crack plane are required; therefore, 2D images of cross sections of concrete specimens are frequently studied; however, the fibre orientation can only be roughly determined from such sections [31,44,45], and the measurement of the embedded length is not possible at all [22]. A more accurate characterisation is obtained by quantitative analysis of micro-computed tomography (µCT) images that allow for a reconstruction of the whole fibre system in 3D [46][47][48].
This work presents a tensile prediction model for UHPFRC-specimens based on a stochastic model for the 3D fibre system and an extensive single-fibre pull-out study. Input parameters for the fibre model are the distribution of fibre orientations, the fibre volume fraction, and the dimensions (length, diameter) of the fibres. The fibre orientation distribution is modelled by a one-parametric distribution family whose parameter β controls the fibre anisotropy, i.e., the scatter of orientations about the preferred direction. The parameter can be estimated either from a sample of fibres observed in a 3D specimen (3D case) or from fibres crossing a planar section of the specimen (2D case).
For fitting the fibre model, a UHPFRC-specimen from a previous study [13] was scanned by micro-computed tomography to reconstruct and analyse its 3D fibre system. This way, embedded length and orientation of the fibres intersecting the crack observed in a mechanical test of the specimen could be determined. A model for the single-fibre contributions to the tensile behaviour is obtained by fitting piecewise linear functions to the force-slip curves observed in single-fibre pull-out tests. Following [34,49], the predicted tensile curve is the sum of all single-fibre pull-out force contributions of crack-crossing fibres. In practice, the tensile behaviour of individually embedded fibres (as used in the pull-out tests) may show deviations from the behaviour of fibres in a larger UHPFRC sample [31]. To take this effect into account, scaling and shifting parameters are introduced, which allow for a calibration of the prediction model to stress-strain curves recorded in experimental tensile tests. Randomising some of the model components yields more realistic curve shapes than a deterministic model.
Our main contribution is the introduction of a stochastic fibre model with a fibre orientation distribution that goes beyond the simple assumptions often considered in the literature (all fibres aligned in one orientation or every fibre orientation is equally probable). Furthermore, we discuss the widespread but flawed assumption that the fibre orientation distribution of crack-crossing fibres coincides with that of the entire 3D fibre system. We present correct and easy parameter estimators for our orientation distribution model. Our fibre model serves as an input to the tensile prediction model, allowing us to simulate the stochastic variability and uncertainty in the tensile prediction. Additionally, certain production parameters such as the fibre volume fraction or the orientation distribution can be varied in the model allowing for a prediction even for cases where no experimental data are available. Apart from the randomness of the fibre geometry, we also randomise the model for single-fibre force contributions. This provides a further means of incorporating uncertainty in our prediction approach and results in more realistic shapes of the stressstrain curves.

Production and Characterisation of the UHPFRC-Specimens
We consider a UHPC-mixture with a maximum grain size of 1 mm. Details on the specification can be found in Table 1 (or see the description of mixture M02 in [13]). Concrete was produced by using the Eirich-Intensive Vacuum mixer of 5 L volume capacity. The same procedure and mix regime as given in [13] were applied. The fresh concrete tests showed similar spread flow, bulk density, and void content values as reported in [13]. For reinforcement, fibres with an ultimate tensile strength of 2800 MPa and elasticity modulus of 200,000 MPa were used. The fibre volume fraction (V V ) was chosen as 2%, and the fibre aspect ratio was l f /d f = 12.5/0.2 [mm/mm]. Specimens of size 40 × 40 × 160 mm 3 were produced. The casting of the specimens was carried out from one side of the formwork, reproducing the configuration M02F2s02 from [13]. The specimens were cured in a climate chamber for 28 days.

Imaging and Image Segmentation
Specimen M02F2s02 from [13] was scanned by micro-computed tomography (µCT) at the Fraunhofer Institut für Techno-und Wirtschaftsmathematik in Kaiserslautern, Germany. To reduce grey value variations in the images, the cubic specimen was placed in a cylindrical UHPC shell during the scanning process. The CT tube was a Feinfocus FXE 225.51 with a maximum acceleration voltage of 225 kV and maximum power of 20W. A Perkin Elmer detector XRD 1621 with 2048 × 2048 pixels was used. The tube voltage was 190 kV, the target electricity 65 µA, and the power 12 W. Tomographic reconstructions were obtained from 800 projections. The specimen of size 40 × 40 × 160 mm 3 corresponds to a reconstructed image of 441 × 441 × 1766 voxels with a voxel edge length of 90.6 µm. The fibre system was segmented from the grey value image as described in [13].
For calculating the tensile force contribution of single fibres in a loaded composite, the individual fibres have to be separated in the segmentation. Due to the coarse image resolution, the space between touching fibres is not sufficiently resolved. Hence, labelling the connected components of the fibre system leads to the formation of fibre clusters. For separating fibres in the clusters, a particle separation based on the watershed transform was applied [50]. During this procedure, fibres were split into segments, which were merged manually to reconstruct the single fibres.
After CT scanning, a four point bending test was performed on the sample, see [13] for details. During this test, the specimen developed an approximately planar crack parallel to the xy plane in the CT image. The location of this crack was identified in the CT scan. Due to the complexity of the single-fibre segmentation, the analysis was restricted to the vicinity of the crack plane such that orientation and embedded length of all fibres crossing the crack could be determined.

Tensile Tests
Three concrete samples were tested experimentally in uniaxial tensile tests. The tensile test regime is as in [15] and briefly summarised in the following. For testing, the specimen was placed in the gripping jaws with a contact area of 40 × 60 mm 2 and fixed by six bolts, which were pulled by a moment of 110 Nm. This pull moment was determined as the maximum such that no cracks occurred in the clamping area. The specimen was fixed in the pull machine by two nuts. In this setup, a field of 40 mm length located in the centre of the specimens was tensioned. The tests were carried out in a displacement-controlled manner. A low load rate of 0.1 mm/min was chosen. The lengthening of the 40 mm field was measured by using two extensometers. The test was stopped as soon as the lengthening of the field reached 2 mm. During the uniaxial tensile tests, load-lengthening curves were recorded.
It is known that eccentricity occurs in the internal resisting forces due to a non-uniform distribution of fibres in the cross-section. Nevertheless, the stress distribution over the cross-section was assumed to remain uniform during the initial cracking phase and after crack localisation. Thus, the uniaxial equivalent stress reads σ = F A , where F is the force value and A is the area of the cross-section. The strain reads ε = ∆L L , where the measured lengthening ∆L is divided by the tensioned length L = 40 mm.
Two strength values were computed for every specimen: the elastic post-crack tensile strength (σ el ), which corresponds to the force at the end of the linear phase in the curve (limit of proportionality) in sense of [51], and the ultimate post-crack tensile strength (σ ult ), which corresponds to the maximum force reached.

Modelling Single-Fibre Pull-Out Curves
For performing single-fibre pull-out tests, individual fibres were embedded in a concrete slab. The same concrete mixture and fibres as described in Section 2.1 were used. The length l e of the fibre part embedded in the concrete and the inclination angle θ of the fibres with respect to the pull-out direction were varied. The embedded length l e was chosen as one half, one third or one sixth of the fibre length, i.e., l e is l f /2 = 6.25 mm, l f /3 ≈ 4.17 mm or l f /6 ≈ 2.08 mm. The inclination angle θ was varied from 0 • to 80 • in steps of 10 • . For testing, the free end of the fibre was clamped between two metal jaws of the testing machine and pulled out while rigidly fixing the concrete slab. During the procedure, force-slip curves reporting the force P applied and the slip s were recorded. At least six fibres were tested for each combination (θ, l e ) resulting in more than 162 single-fibre pull-out curves. Figure 1 shows the force-slip curves P i (s, θ, l e ), i = 1, . . . , 6, of the six tested fibres with θ = 10 • and l e = l f /2. We denote byṖ(s, θ, l e ) the median curve, which is the point-wise median of the six single-fibre pull-out curves for a fixed combination (θ, l e ), i.e., P(s, θ, l e ) = med(P 1 (s, θ, l e ), . . . , P 6 (s, θ, l e )), s ∈ [0, l e ]. (1) In the next step, we fit a piecewise linear model to the observed force-slip curves. In the literature, several approaches (bilinear, trilinear, exponentially decreasing) for modelling single-fibre pull-out curves have been suggested, see [31]. Based on the results of our single-fibre pull-out tests, we decided to use a three-phase model. Phase I represents the linear elastic part of the curve up to the yield strength P el that is reached at slip s el . Phase II is the nonlinear part up to the ultimate force P ult at slip s ult . For simplicity, this part is also described by a linear model. Phase III is a linear descending branch up to complete fibre pull-out at slip s tot , i.e., the last recorded slip value. See Figure 2 for an illustration of the model.  For fitting the model, the force and slip values derived from the six fibre-pull-out tests per combination (θ, l e ) are averaged. We use the median values P el,med (θ, l e ) = med(P el,1 (θ, l e ), . . . , P el,6 (θ, l e )) P ult,med (θ, l e ) = med(P ult,1 (θ, l e ), . . . , P ult,6 (θ, l e )) s el,med (θ, l e ) = med(s el,1 (θ, l e ), . . . , s el,6 (θ, l e )) s ult,med (θ, l e ) = med(s ult,1 (θ, l e ), . . . , s ult,6 (θ, l e )) s tot,med (θ, l e ) = med(s tot,1 (θ, l e ), . . . , s tot,6 (θ, l e )).
The numerical values are given in Appendix A, Table A1. When performing a robust two-factor ANOVA, equality of means is rejected for each of the five characteristics. One exception is P ult where the hypothesis of equal means over angles cannot be rejected at the 5% level (p-value = 0.083). In spite of this finding, we use individual values for all (θ, l e ) groups for all five characteristics rather than merging groups for P ult .
In summary, the model reads Formulas for and estimates of the parameters p 1 , p 2 > 0, p 3 < 0, and r 1 , r 2 ∈ R are summarised in Appendix A, Table A2. Figure 3 illustrates the model for θ = 10 • and l e = l f /2, l f /3, l f /6.

Prediction Model for Tensile Stress
Our model is based on the following assumptions (see also [16,49]): • The UHPC matrix is a statistically homogeneous material [52], (p. 28).
• Fibres are straight with cylindrical geometry. The fibres have fixed length l f and diameter d f . • The spatial distribution of the fibre positions in the UHPC is statistically homogeneous [52], (p. 28). • The fibre orientation is random following a given probability density function (p.d.f.) p.

•
The specimen is uniaxially strained. • During loading, the specimen develops a planar crack C of width w and area A C orthogonal to the tension axis, see Figure 4. • With growing crack width, the shorter fibre end is pulled out of the concrete matrix. The longer end is not affected. • The fibres behave linearly elastic. • The matrix deformation and the Poisson effect of the fibres during pull-out are neglected. The fibre-matrix bond is frictional.
If the strain ε exceeds ε r , the strain at the uniaxial tensile strength of the UHPC, a crack starts to form. The crack width w is given by with L = 40 mm and ε r = 0.00087. Fibres crossing the crack are divided into the two segments to the left and to the right of the crack plane. We denote the length of the shorter segment by l e , see Figure 4, such that l e ∈ [0, l f /2]. With increasing crack width w, the shorter fibre end is pulled out of the concrete until the fibre becomes detached at w > l e . . UHPC matrix with one fibre intersecting a planar crack C (assumed to be embedded in the (X,Y)-plane) at crack width w = 0. The crack divides the fibre in two parts. The embedded length l e is the length of the shorter (green) fibre part. The angle between fibre and tension axis (assumed to be the Z-axis) is denoted by θ.
For predicting the stress evolution with increasing crack width w, the single-fibre pull-out forces P(s, θ, l e ) of the fibres crossing, the crack have to be taken into account. Due to the assumption that only the shorter fibre end is pulled out of the concrete, the slip s corresponds to the crack width w. Furthermore, l e in the single-fibre pull-out tests corresponds to the embedded length l e . Thus, we use s = w and l e = l e in (2).
Assume that the planar crack C is crossed by N fibres with given inclination angles and embedded lengths (θ 1 , l e,1 ), . . . , (θ N , l e,N ). As experimental single-fibre pull-out curves are only available for a discrete set of θ and l e values, the forces P(w, θ, l e ) of fibres in C are binned into classes as given in Table 2. Table 2. Binning of embedded length l e and inclination angle θ.

Class
Interval P(w, θ, l e ) Fibres crossing the crack counteract the opening of the crack. According to [34,49], the composite stress (or mean resistance force per unit area) at crack width w is obtained by where p c (θ, l e ) denotes the joint probability density of inclination angle and embedded length of crack-crossing fibres. λ c is the mean number of fibres per unit area in C.
We approximate σ ct (w) by first replacing the integral in (4) by a sum over all fibres observed in the crack. In the second step, angles and embedded lengths are binned such that model (2) can be applied.
Here, N T i ,L j denotes the number of fibres in C with θ k ∈ T i and l e,k ∈ L j , k ∈ {1, . . . , N}, i ∈ {1, . . . , 9}, j ∈ {1, 2, 3}.P is a prediction of P by the median curve or by the trilinear model. We denote the prediction based onP =Ṗ by σ med ct (w) and the prediction based oñ P = P tri by σ tri ct (w).

Stochastic Fibre Model
In the following, we introduce a stochastic fibre model to generate 3D fibre systems of virtual specimens. Their tensile behaviour can then be predicted with our stress prediction model.
The fibre system is modelled by a Boolean model [53] as follows: The positions of fibres are indicated by their midpoints, which are modelled by a Poisson point process. That is, the number N of fibre midpoints (x, y, z) in a given volume V follows a Poisson distribution with parameter λ · V, where λ > 0 is the mean number of fibres per unit volume. Locations of fibre midpoints are drawn independently from a uniform distribution on the volume of interest.
The orientation of a fibre is independent of its location and can be described in spherical coordinates with co-latitude angle θ ∈ [0, π/2) and longitude angle ϕ ∈ [0, 2π). Assuming the Z-axis to be the tension axis, θ corresponds to the inclination angle of the fibre with respect to the tension axis. The longitude angle ϕ does not influence the force contribution of a fibre. Hence, we consider ϕ to be uniformly distributed on [0, 2π) such that the probability density function of the fibre orientation is a function p(θ) depending only on θ.
We assume that θ follows a one-parametric orientation distribution described by the p.d.f see [50,[54][55][56] for details and applications. The anisotropy parameter β controls the alignment of the fibres. For β = 1, the fibres are isotropically oriented. For decreasing β, the fibres tend to be aligned along the Z-axis, see Figure 5. For β > 1, the fibre orientations are concentrated in a plane. This case is not considered here.
(a) (b) (c) For fitting the model to a given concrete sample, the characteristics λ and β can be obtained from µCT images, see [50,55,57]. Due to the low fibre volume fraction, we assume that overlap of fibres in the model is negligible. Hence, the fibre intensity λ can be computed from the fibre volume fraction V V , the fibre length l f and cross-sectional area A f via If a single-fibre segmentation is available, an estimate of the anisotropy parameter β can be determined from the sample θ 1 , . . . , θ N of inclination angles by using the maximum likelihood method as described in [54]. An estimate of β based on the method of moments is given by (see Appendix B)β The Boolean model described so far yields a stochastic fibre system in 3D. The prediction model outlined in Section 2.5 requires only inclination angles θ and embedded lengths l e of the fibres intersecting the crack plane C. Due to the spatial homogeneity of the model, C can be assumed to be contained in the (X,Y)-plane. In the Boolean model, fibre position and orientation are independent. Hence, l e and θ are independent. Thus, the joint p.d.f. of θ and l e for crack-crossing fibres fulfils p c (θ, l) = p c (θ)p l e (l) where p c (θ) is the p.d.f. of inclination angles θ of fibres crossing the crack and p l e (l) is the p.d.f. of l e . The spatial homogeneity of the Boolean model implies that l e is uniformly distributed on [0; l f /2]. The planar characteristics are related to the spatial characteristics as follows [49,50,53] p l e (l) = 2 l f , l ∈ [0; l f /2] (10) where λ c is the expected number of fibres per unit area in C.
Note that the difference between (7) and (11) can be explained as follows [49,58]: The p.d.f. p(θ) takes every fibre in the reinforced concrete into account whereas p c (θ) accounts only for fibres intersecting the crack plane. A fibre is more likely to hit the crack plane if its inclination angle θ is small. The assumption that the distribution p c (θ) is equal to p(θ) is a common misunderstanding that leads to incorrect interpretations of measurement results (see [22,31,59]).
If a single-fibre segmentation of fibres in a crack is available, an estimate of the anisotropy parameter β can be determined from the sample θ c,1 , . . . , θ c,N of inclination angles of crack-crossing fibres by using a modified version of the maximum likelihood estimator given in [54]. An estimator of β based on the method of moments is given bŷ see Appendix B for details. The fibre orientation coefficient η θ is a widely established characteristic for the distribution of fibre orientation in a crack [22,31,33,50,[60][61][62]. η θ is the second moment of the angular deviation of the crack-crossing fibre from the tension axis. Here, it is possible to calculate η θ in closed form: For an isotropic fibre orientation (β = 1), which is assumed in several studies [19,22,31,34], it follows that λ c = 1 2 λl f , p c (θ) = 2 sin(θ) cos(θ) = sin(2θ) and η θ = 1 2 . In the extreme case that all fibres are aligned along the tension axis (β = 0), it follows that λ c = λl f , and η θ = 1.

Image Analysis
The crack observed in the imaged sample is approximated by a plane whose position is identified in the CT scan, see Figure 6a. Among the segmented fibres, N = 598 cross the crack plane. Inclination angles of the individual fibres are determined by using partial second derivatives as described in our former study [13]. Fibre lengths can be measured by the maximal Feret diameter [63]. Note that the Feret diameter of a particle is defined as the minimal distance between two parallel planes that are orthogonal to a specified direction and enclose the particle. The maximal Feret diameter is obtained by maximization over all directions. By intersection with the crack plane, each fibre is split into two segments. The embedded length l e is the maximal Feret diameter of the shorter one of these segments. From the inclination angles θ c,i , we obtainβ c,mom = 0.22.
By Equation (12) with β = 0.22 the expected number of fibres in the cross section is A C · λ c = 834. The observed number N = 598 is well below that value. There are several possible explanations. In [13], the fibre content at the crack location was observed to be about 5% lower than the theoretical value of V V = 0.02 (see [13] Figure 12g at the approximate crack position in Slice 641). Additionally, according to the master datasheet of the fibres used, both, the diameter d f and the fibre length l f can deviate by 10%. In Figure 7, the distributions of observed inclination angles and embedded lengths are compared with the fitted p.d.f.s p c (θ) (with β = 0.22) and p l e (l). The inclination angles are well fitted. A chi-square goodness-of-fit test does not reject the hypothesis that the inclination angles follow a distribution with p.d.f. p c (θ) with β = 0.22 (p-value = 0.1244). The distribution of the embedded length has an unexpected peak around 4.5 mm. The hypothesis that the embedded length are uniformly distributed is rejected by a Kolmogorov-Smirnov goodness-of-fit test (p-value < 10 −4 ). A possible explanation is that short fibre segments (whose diameter in the µCT image corresponds to only 2 to 3 voxels) were missed. It should be noted that we are not aware of any practical way of estimating l e from a 2D slice as was also mentioned in [22]. A fibre system simulated as a realisation of the Boolean model with β = 0.22, l f = 12.5 mm, d f = 0.2 mm, and V V = 0.02 is shown in Figure 6b.

Prediction of Tensile Stress
To predict the composite stress σ ct using Equation (6), we need N T i ,L j , the number of crack-crossing fibres with (θ, l e ) in class (T i , L j ), i = 1, . . . , 9, j = 1, 2, 3. Based on the single-fibre segmentation, we derive N T i ,L j as given in Table 3.
Experimental tensile tests were carried out on three specimens as described in Section 2.3. The ultimate tensile stress and the corresponding strain of the experimental tensile curves are given in Table 4. The observed tensile curves are used to calibrate our prediction.

Prediction of Tensile Stress Based on Median Curve
Here, we use the median curve to predict σ ct by σ med ct . Figure 8a shows the prediction as well as the three experimental tensile curves. We see that the prediction overestimates the stress. The predicted ultimate tensile stress is 18.56 MPa which is reached at strain ε med ult = 0.0230 [mm/mm]. To calibrate the curves, we follow [19,22,31,64] and introduce suitable rescaling factors to match characteristic stress-strain points (i.e., ultimate tensile stress σ ct,ult or yield stress σ ct,el ) of experiments and predictions.
In the first step, we scale σ med ct by a stress scaling factor S σ ult such that the predicted and mean experimental ultimate stress coincide, i.e., S σ ult σ med ct,ult = σ ult,exp . This is achieved for S σ ult = 0.61, see Figure 8b. In the next step, the locations of the maxima have to be matched. To this end, we additionally scale the strain axis by a strain scaling factor S ε ult . The choice S ε ult = 0.42 yields S ε ult ε med ult = ε ult,exp such that the maxima are closer together (see Figure 8c); however, the slope of the curve after the maximum does not fit the experimental curves. A remedy is to restrict scaling of the strain axis to the region prior to the maximum, that is, where ε med tot is the last available strain value in the median curve. This way, both, the maximum location and the stress at the final point ε med tot are matched. The result shown in Figure 8d indicates that the slope after the ultimate tensile stress is reduced due to the fixation of ε med tot .
The new parameters (p • 1 , p • 2 , p • 3 , r • 1 , r • 2 ) are calculated as in Appendix A2. The prediction is given in Figure 9. The factors are chosen such that σ tri ct,ult = σ ult,exp and ε tri ct,ult = ε ult,exp . In the interval between yield and ultimate stress, the predicted curve is straighter than the experimental stress curves, which can be explained by the use of a linear model. Furthermore, the maximal peak is more pronounced than in the experimental tests. This may be due to the replacement of the individual single-fibre pull-out forces P i by average forces P tri obtained from the trilinear model. To overcome this problem, we propose a model based on randomised single-fibre contributions in the following section.

Prediction of Tensile Curves Based on Randomised Trilinear Model
In order to model the random variations in the tensile curves, we consider s ult,med as random following a uniform distribution and add normally distributed residuals to the yield and the ultimate force. This way, individual functionsP are simulated for each fibre in the crack. The range of values of the required random variables is inferred from the single-fibre pull-out experiments.
M = 10 predictions based on the fibre system observed in the CT image are shown in Figure 10a. The predictions reproduce the stress profile better than the deterministic predictions.
Prediction study part 2: Stress prediction of virtual specimens with varying production parameters We use the stochastic fibre model from Section 2.6 to generate fibre systems of virtual specimens. The model parameters are determined from the crack-crossing fibres of the scanned specimen. From each model realization, we derive the values of (θ, l e ) for all fibres intersecting a virtual crack plane C. Figure 10b shows a tensile prediction of ten virtual fibre systems when using the randomized prediction scheme described above.
Predictions for virtual samples with modified orientation distribution and volume fraction V V are shown in Figure 10. Note that the fibre parameters d f and l f were not varied since the single-fibre pull-out tests were restricted to fibres with d f = 0.2 mm and l f = 12.5 mm. Figure 10c shows that a fibre orientation along the tensile axis increases the tensile stress compared to an isotropic fibre orientation. This is in accordance with [22,31].
In particular, this plot shows that the common assumptions of completely aligned or isotropic fibres result in quite different predictions. Figure 10d reveals that the influence of the volume fraction also behaves as expected: Fewer fibres decrease and more fibres increase the tensile stress. Figure 10. Prediction of M = 10 tensile curves (a) based on the CT scanned fibre system summarised in Table 3

Conclusions
In this study, we present a prediction model for tensile behaviour of UHPFRCspecimens based on statistical information on the fibre system and extensive single-fibre pull-out tests. The model is calibrated by comparing the prediction to the results of experimental uniaxial tensile tests. We introduced a stochastic fibre model that allows for the generation of fibre systems of virtual specimens, which can be used for the prediction of the tensile behaviour. Through experimental and theoretical investigations, the following conclusions are drawn: • The fibre system in the concrete is modelled by a Boolean model of straight cylinders. The fibre orientation distribution is represented by a one-parametric distribution family. Its parameter β controls the anisotropy of the fibre orientation. Both, the case of total alignment and isotropy are included in the model as special cases. The widely used reference value η θ (fibre orientation factor) can be calculated directly from β.
In particular, the relationship of fibre orientation of the full 3D fibre system and of crack-crossing fibres is examined. • We emphasised the difference between the p.d.f. for the fibre orientation of all fibres in the reinforced concrete and the p.d.f. for the fibre orientation of fibres intersecting a crack plane. The assumption that the p.d.f.s are equal is a common misunderstanding that leads to incorrect interpretations of measurement results (see [22,31,59]). • Estimators for β are presented and analysed. We recommend the method of moments estimatorβ c,mom . It estimates accurately, no numerical method is required (such as a gradient descent algorithm in [55]) and only the fibre angles in 2D cross-sections are needed. • In contrast to 2D imaging methods, µCT allows for a determination of the embedded length of fibres in cracks since the whole 3D fibre system is observable. In particular, this allows for the validation of the fundamental assumption of a uniformly distributed embedded length l e for the widespread tension-softening model in [16]. Furthermore, using µCT overcomes the accuracy disadvantage of fibre orientation estimation in 2D sections as discussed in [22]. • The presented prediction model uses the stochastic fibre model combined with a statistical analysis and stochastic modelling of the single-fibre pull-out tests. Predictions are calibrated by using experimental tensile curves. Running the prediction on virtual specimens with varied production parameters led to reasonable stress-strain curves.
More research with additional samples is needed to validate the model and investigate its robustness to changes in the production parameters.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: UHPFRC ultra-high performance fibre-reinforced concrete µCT micro-computed tomography p.d.f.
probability density function area of the specimen cross-section C planar crack A C area of crack C w crack width s slip θ inclination angle l e , l e embedded length in single fibre pull-out tests, in cracked UHPFRC s el (θ, l e ) slip value of P el (θ, l e ) s ult (θ, l e ) slip value of P ult (θ, l e ) s tot (θ, l e ) last slip value s el,med (θ, l e ) median of s el (θ, l e ) s ult,med (θ, l e ) median of s ult (θ, l e ) s tot,med (θ, l e ) median of s tot (θ, l e ) P(s, θ, l e ) single-fibre pull-out force P el (θ, l e ) tensile force at the end of the linear phase P ult (θ, l e ) ultimate forcė P points-wise median curve P tri trilinear curve P el,med (θ, l e ) median of the six fibre pull-out forces P el (θ, l e ) P ult,med (θ, l e ) median of P ult (θ, l e ) p i , r j parameters of P tri , i = 1, 2, 3, j = 1, 2 T i , L j classes for θ, l e , i = 1, . . . , 9, j = 1, 2, 3 normal distribution with mean a and variance b 2 I ult (θ, l e ) interval giving the range of s ult,med (θ, l e ) over all θ but a fixed l e s * ult,med (θ, l e ) random value based on interval I ult (θ, l e ) sd P ult,med (θ, l e ) sample standard deviations of P ult,med (θ, l e ) sd P el,med (θ, l e ) sample standard deviations of P el,med (θ, l e ) r * P ult,med (θ, l e ) random values based on sd P ult,med (θ, l e ) r * P el,med (θ, l e ) random values based on sd P el,med (θ, l e ) Appendix A Table A1. Median values s el,med (θ, l e ), P el,med (θ, l e ) (standard deviations given in brackets), s ult,med (θ, l e ) (minimum and maximum values given in brackets), P ult,med (θ, l e ) (standard deviations given in brackets), s tot,med (θ, l e ). No values for sd P el,med and sd P ult,med were determined for the combination (θ, l e ) = (80 • , l f /6), as only one single-fibre pull-out curve could be observed.  Table A2. Parameter values for P tri (w, θ, l e ). No value for p 1 was determined for the combination (θ, l e ) = (80 • , l f /6), as no phase I could be observed.

Appendix B
We consider a fibre system as a realisation of a set of random line segments. The random angular deviation between a line segment and the Z-axis is described by its inclination angle θ. A random inclination angle in space is denoted by Θ and a random inclination angle of a fibre intersecting the (X,Y)-plane by Θ c . Θ follows a distribution with density p(θ) = β sin(θ) (1 + (β 2 − 1) cos 2 (θ)) 3 2 , θ ∈ 0, π 2 , β > 0.
The difference between p(θ) and p c (θ) arises from the fact that p c (θ) only considers the fibres intersecting the (X,Y)-plane, whereas p(θ) considers every fibre of the fibre system [49]. In a realisation with N line segments, we denote a sample of spatial inclination angles by θ 1 , . . . , θ N and a sample of inclination angles of fibres intersecting the (X,Y)-plane by θ c,1 , . . . , θ c,M , M ≤ N.

Method of Moments Estimatorsβ mom ,β c,mom for β
The method of moments (mom) is based on the approximation of sample moments with theoretical moments (via the law of large numbers (LLN)). If the theoretical moments are functions of a parameter, the sample moments are functions of the parameter's estimate. Here, the parameter of interest is the anisotropy parameter β. Assume that Θ 1 , . . . , Θ N i.i.d. with density p such that for suitable function g we have E(Θ k i ) = g k (β) < ∞. Then an estimatorβ mom can be obtained as solution of Θ k N = 1 N ∑ N i=1 Θ k i = g k (β mom ). We apply the cosine to the angles in the following. Thus, cos(Θ 1 ), . . . , cos(Θ N ) are i.i.d. with E(cos(Θ i )) = π 2 0 cos(θ)p(θ) dθ = 1 1 + β (A8) and g(β) = 1 1+β , β > 0. Therefore,β