Finite element modelling of UHPC under pulsating load using X-ray computed tomography based fiber distributions

The benefits of including fibers in ultra-high performance concrete (UHPC) are attributed to their good bond with the matrix and, hence, an optimal utilization of their properties. At the same time, though, fiber reinforcement may contribute to anisotropy in the composite material and induce weak areas. The influence of the fibers’ orientation on the material properties is a matter of current scientific discourse and it is known to play a vital role in structural design. In the case studies presented herein, mechanical laboratory tests using pulsating load regimes on UHPC with a strength of more than 200 MPa were simulated by use of finite element models. The orientations of the fibers were measured for each test sample prior to failure using an X-ray computed tomography (CT) scanner, and these orientations are explicitly implemented into the model. The paper discusses the methodology of merging data retrieved by CT image processing and state-of-the-art FE simulation techniques Moreover, the CT scanning was carried out throughout the testing procedure, which further enables the comparison of the mechanical tests and the FE models in terms of damage propagation and failure patterns. The results indicate that the overall fiber configuration and behavior of the samples can be realistically modelled and validated by the proposed CT-FE coupling, which can enhance the structural analysis and design process of elements produced with steel fiber reinforced and UHPC materials.


Introduction
With concrete being currently the most consumed material by humanity after water [1], but also one with a significant carbon footprint and one of the most demanding in terms of energy resources [2,3], it is becoming obvious that an increase in its performance is vital. A significant improvement is the production of concrete materials with very high compressive strength (e.g. higher than 120 MPa), which is referred to as ultra-high performance concrete (UHPC) or reactive powder concrete and finds application where challenging construction material characteristics are required, for example in high-rise, industrial, and infrastructure projects. This type of compressive strength can be achieved through the use of generally finer aggregates and lower water/cement ratios than normal concrete, and the addition of appropriate cementitious and additive constituents [4]. The addition of dispersed short steel fibers in UHPC -but also in normal concrete -is a commonplace engineering decision (hereafter UHPC is assumed de facto to be steel fiber reinforced), since it leads to constrained concrete crack formation and propagation, which is then associated with various design and performance advancements: ductility and toughness, freeze-thaw and shrinkage resistance, sustained load and flexural creep resistance, as well as improved dynamic fatigue and impact response [5,6].
However, the addition of fiber reinforcement can have critical disadvantages and even hinder the performance of concrete, since it can contribute to the anisotropy of the mixture and induce weak areas. The fiber distribution strongly depends on the casting conditions of the fresh mix as well as geometric constraints, and it has been shown that different casting methods for the same mix, and hence different prevailing fiber orientations, can lead to significant differences in structural performance [7][8][9]. The review by [10] of a broad range of research attempts to predict the fiber distribution and orientation confirms that this has been possible in individual areas of applicability, but there is no unanimously accepted model capable of considering all the possible variations such as type of fiber, fiber geometry, concrete rheology, fiber stiffness, fiber dosage and the compacting method, but as a general consensus, the alignment of fibers in a flowing medium tend to adopt the flow lines. Some concepts have also been proposed regarding control of the fiber alignment in fresh concrete by use of mechanical means [11], flow techniques [12], as well as magnetic fields [13,14].
Beyond the attempts to control the fiber alignment quality, a variety of methods is available in order to implicitly evaluate the distribution and orientation of steel fiber reinforcement in hardened concrete. A straightforward method is by observation and recording the fibers visually in a fractured or cut plane of a sample [15]. This method can be automatized by use of image analysis tools and software, which can provide an estimate of the number and orientation of fibers in cut fiber reinforced concrete specimens by analyzing the cross section shapes of the fibers [16,17]. A range of further non-intrusive methods are proposed by various researchers. Inductive methods make use of the magnetic resistivity of steel fibers and concrete by subjecting a concrete specimen to an electromagnetic field and measuring the magnetic flow [18,19]. The attribute of electrical resistivity of the fibers has been utilized in order to assess the fibers' orientation, their density, and their degree of segregation, e.g. by use of impedance measurements [20][21][22]. The possibility of microwave permittivity measurements of steel fiber reinforced elements has been shown in [23].
Unlike the methods mentioned above, which implicitly address the issue of fiber orientation, X-ray computed tomography (CT) can be used to achieve an imaging and full record of the fibers in their exact location and distribution within a cast or cut concrete specimen. X-ray computed tomography has been applied to concrete research problems for several decades [24], and it is particularly well suited for investigations on fiber reinforced material, especially when combined with measured mechanical response data that can be related to individual fiber spatial distribution and directional orientation [25]. Especially for applications with UHPC, [25] demonstrates a procedure to record the orientations of fibers in the matrix and transform them to a probability distribution by means of CT. Previous research has also achieved CT scanning at different deformation increments during unconfined compression, cylinder splitting, and triaxial tests [26][27][28][29][30], which facilitated the measurement of cracking initiation and propagation in a single specimen as a function of the load and deformation state.
Simultaneously, non-linear structural modelling of fiber reinforced concrete is vital for the efficient design of structures, since it allows for realistic prediction of the damage and fracture behavior of the construction material. Several attempts have been made to model the performance of fiber reinforced concrete and there is a plethora of tools available to predict the responses of such materials under various loading conditions. Most approaches are based on the cohesive crack model [31] which makes use of a r -w or stress-CMOD (Crack Mouth Opening Displacement) response of the material at the mesoscale. In this case, the fiber orientation and distribution are not taken directly into consideration and the r -w response is obtained experimentally based on mechanical testing on specimens. Despite these simplifications, the cohesive crack model currently remains the standard way to characterize materials and design concrete structures in practice [32].
In finite element (FE) simulations, the varying effect of the fiber orientations has been accounted for by non-isotropic material laws [33] and spatial variations of the fiber volume fraction [34]. An approach of explicitly modelling fibers in cementitious materials based on discrete inclusion of fibers and to account for the contribution of each single fiber based on its location and alignment has been proposed in [35], and since then extensively studied. In [36], a method for replacing the explicit fibers with discrete forces integrated with a background mesh was proposed and validated. In [37], a two-phase model, where fibers and the concrete are modelled separately, was proposed. The individual fibers' load-displacement performance was evaluated based on fiber pull-out tests and simulated by a multi-linear, axial load-displacement law, while the concrete matrix performance was simulated with a smeared crack constitutive law. The fiber dispersion was then generated randomly by a Monte Carlo algorithm with theoretical distribution functions without test-based information [38], but the FE analyses compared well with actual monotonic quasi-static material tests. Using this modelling philosophy, [39] shows the sensitivity of the simulations to the variations of the fiber orientations and successfully validates an envelope of multiple model results against an envelope of bending tests with normal strength fiber reinforced concrete. A similar approach is used in [40], where prisms are modelled using nonelastic constitutive laws for the fiber and the concrete matrix, and a cohesive interface function is validated against beam bending and direct tensile tests for concrete with a compressive strength of up to 110 MPa with random fiber distributions. A mesoscopic modelling approach for cement composites, combining models for discrete cracks using the strong discontinuity approach, explicit fibers as axial truss elements and individualized fiber bond spring-models is also proposed in [41]. A further mesoscopic FE simulation technique with explicit description of the concrete, randomly distributed fibers, and a fibermatrix bond is demonstrated in [42].
A multi-level model capturing the single fiber, the crack plane, and the continuum/component has been proposed by [43]. This study accounted, at the first level, for an analytical pull-out model and a castinginduced anisotropic fiber distribution. At the second level, it accounted for a stress transfer across a crack (through integrating the crack-bridging contribution of randomly aligned fibers) and material cohesion (dependent on the tensile strength of plain concrete). The analysis at continuum level provided very good agreement not only with quasi-static tests at the level of material and simple components testing [43], but also with preliminary results of non-reinforced cyclic 3-point bending tests as part of an ongoing model validation for cyclic/dynamic testing [44].
As dictated by the current research background discussed above, CT methods can accurately record steel fibers in concrete materials and state-of-the-art FE simulation techniques allow for a discrete representation of the fibers and their mechanical contribution in non-linear models of small concrete components. To the knowledge of the authors, a merging of these two possibilities has not been presented in the past, and it is attempted here for the first time in international literature, in terms of a case study. A further novelty of this investigation is that focuses on FE modelling of UHPC with a strength of more than 200 MPa, and on behavior under pulsating loading. The study demonstrates a methodology for the FE modelling of UHPC components by implementing actual fiber orientation information based on CT measurements and image analyses of actual samples. The mechanical tests and FE simulations follow the same pulsating loading regime, where the simulated tests are compressive and double-punch tests. The CT imaging takes place incrementally, before and after each load step and it, hence, captures not only the exact location and orientation of the fibers, but also the material damage propagation. The discussed FE modelling is then validated in terms of fiber configuration (geometrical validation) and in terms of compressive and double-punch pulsating load-displacement response and failure pattern (mechanical validation). The described methodology of coupling CT scanning with FE modelling can then enhance the structural analysis and design process for large-scale structures and structural components made of UHPC.
Following this introduction, the UHPC material ''Cor-Tuf'' that is implemented in the present study is described. Then, the experimental investigations in combination with CT scans are discussed. The modelling procedure, including the techniques of incorporating the fibers discretely in the model and carrying out the FE simulations, is presented in detail. The results and a comparative assessment of the specimens' mechanical response in the experimental and numerical investigations lead to the main conclusions of this paper.

Material
The material used in study is designated as Cor-Tuf, which is the name of a family of UHPC mixes developed at the U.S. Army Engineer Research and Development Center (ERDC). Cor-Tuf is distinguished by a high compressive strength, which generally ranges from around 190 MPa to 250 MPa and with mixture proportions depending on desired fresh and hardened properties required for individual projects [45]. The version of Cor-Tuf used herein is the standard mix design (also known as Cor-Tuf1), because it has had a good track record of performance and has been extensively characterized [45,46]. The material has a maximum particle size of approximately 0.6 mm, which is governed by the maximum size of the foundry-grade Ottawa sand used in the mix. The mixture proportions for Cor-Tuf are given in Table 1, and further details of the steel fibers used are provided in [47]. The fibers were hook-ended and had a length of 30 mm and a diameter of 0.55 mm. Table 2 provides representative mechanical values. All values are based on studies reported in [45], except for values resulting from four-point bending or splitting tensile tests, which were taken from [46]. In should be noted that although the direct tension tests and flexural/ splitting tests are different, the resulting values agree well with the general assumptions for the derivation of the direct tensile [48] and the flexural/splitting strength of concrete from its compressive strength [49] for typical normal and high-strength concrete. Moreover, [46] observes that during flexure tests of Cor-Tuf, little or no fiber rupture or straightening of the hooked fiber ends tended to occur during the failure process, which indicates that the matrix material of Cor-Tuf tends to fail around the fibers, leaving the hooked ends intact. Investigations comprising pull-out tests of fibers with the same strength and hook-ended anchorage (3D DramixÓ) but with much larger cross sectional areas have been carried out for high-strength concrete [50] and UHPC [51] with compressive strengths of 68 and 158 MPa, respectively. During these pull-out tests, fibers in both materials showed a strong tendency for the failure to occur as fiber rupture and hence full utilization of the fiber strength.

Experimental set-up
A series of factors related to the X-ray CT system and the mechanical testing apparatus available for use during these experiments dictated the specimen size. Given the high strength of the material, the load limit of the testing equipment for compression tests is reached with cylinders of approximately 70 mm in diameter. The choice of this specimen size also considered the requirement for incremental loading up to failure with CT scans collected after each loading step to capture damage propagation. Furthermore, the quality of the CT images was found to significantly deteriorate for specimens with diameters exceeding approximately 76 mm. Due to these limitations, the nominal specimen diameter was selected to be 70 mm. The cylindrical specimens of Cor-Tuf were cut from cores drilled out of the middle of a thick concrete block in order to avoid constraints on fiber orientation caused by the boundaries of small casting molds and to obtain a concrete performance as closely as possible representative of actual structural members. The original block had a size of 914 by 914 by 556 mm and was cast directly from a concrete mixing truck. Two cylindrical test specimens were used in this investigation, namely one for unconfined compression and one for double punch testing. The specimens used in this study are marked as UC-EXPfor uniaxial compression and DP-EXP for the double punch test and they had the following dimensions: • UC-EXP diameter and height: 70.3 mm and 141.3 mm, respectively • DP-EXP diameter and height: 70.4 mm and 70.5 mm, respectively The unconfined compression (UC) tests were carried out using geometric characteristics and testing procedures similar to those recommended by the procedures in [52]. The double punch (DP) testing method (also referred to as the Barcelona test) is an option for obtaining tension characteristics of concrete materials. Although much less common than direct tension, cylinder splitting, or flexural beam tests, it offers a number of significant benefits related to the demands of both X-ray CT scanning and incremental load application. The test geometry followed the methodology provided in [53]: the concrete cylinder was placed between two cylindrical steel punches centered on its top and bottom surfaces (Fig. 1). Each punch had a diameter equal to the concrete cylinder diameter and the two punches were compressed towards one another during the experiments. This is intended to create a failure mechanism that includes a shear cone under each steel punch and a series of radial cracks that protrude from the center of the concrete specimen out to its circumference. The tests were displacement controlled, conducted using a 979-kN capacity MTS testing machine. Load was recorded using the load cell incorporated into the MTS testing machine and displacement was recorded using external linear variable differential transducers (LVDTs). A photograph of this testing configuration can be seen in Fig. 2. The recording rate for all data was 20 Hz. The load history of both tests included a series of progressively increasing pre-damaging loading cycles, which were only halted once the effective failure of the specimen had occurred. For the DP specimen, additional loading cycles following peak load were also completed to investigate its post-peak ductile behavior. The loading rate was representative of quasi-static loading situations and after most loading cycles, the unloaded specimens were scanned using X-ray CT. For the X-ray CT scanning, the specimens were placed on a rotating table between an X-ray source and an X-ray detector, projecting an X-ray attenuation image of the specimen upon the detector. Once the projected images were recorded during the 360 rotation of the specimen, mathematically-based reconstruction algorithms were used to produce a threedimensional representation of X-ray attenuation within the specimen. The recorded material density, and hence the individual component materials and voids in the specimens, could then be correlated with the recorded X-ray attenuation. This led to the identification of the specimen material structures and a separation and analysis of the individual components with the same order of density (voids, matrix, steel fibers). The scanning system that was used during all X-ray CT experiments was located at the University of Florida. A 225-kV microfocus X-ray tube manufactured by Comet GmbH was used to generate the X-ray beam. The flat panel detector was manufactured by Thales Electron Devices and had an active image sensor area of 285 by 406 mm with an array of 2240 by 3200 pixels and a pixel size of 127 lm. The 3D-CT images were reconstructed using the program efX CT   [58] and all subsequent image analysis was completed in a self-programmed code in MATLAB.

Results of CT scans and mechanical tests
An X-ray CT-based rendering of the fiber characteristics in the immediate vicinity (± 25 image slices) of a series of cross sections are provided in Fig. 3 for each of the specimens along with the definition of the corresponding cross section locations. Figure 4 provides a series of X-ray CT image cross sections showing the internal crack structures following the final loading cycle for each specimen. The entire set of load-displacement curves for each specimen will be shown in Fig. 10, in combination with the numerical simulation results.
Using the X-ray CT data, it was possible to extract detailed information about fiber orientation characteristics using a Hessian-based approach [25], which analyzes fiber orientation at each fiber voxel (i.e. 3D pixel) point. By combining each of these orientation measurements, it is possible to obtain an estimate of the global distribution of fiber orientation within each specimen. The results of these fiber orientation analyses for each sample will be presented in Sect. 4.1.3 along with a comparison to the fiber orientations generated in the numerical models. These results will be presented using a spherical coordinate system, in which orientations are characterized by angles h and U (see Fig. 5left). The angle h represents the azimuthal angle in the x-y plane from the x-axis (in this context, the cylindrical axis of the sample is denoted as the z-axis), with -180°\ h \ 180°. The angle U represents the polar angle from the positive z-axis. Note that since the fibers are assumed to be symmetric about their lengths, a corresponding symmetry condition can be imposed on U, with 0 \ U \ 90°. For more details about the fiber orientation analysis, see [25].
The failure of the specimen under UC was characterized by many small vertical cracks and columnar   formations through both ends (see Fig. 3). This is largely coincident with the anticipated failure pattern of well-performing homogeneous fiber reinforced concrete, since these many small cracks correspond to high energy absorption, as compared to the failure patterns of unreinforced concrete, which are typically characterized by large diagonal fractures. However, the fiber orientation range was not entirely uniform in the specimen, which partially explains some relatively wide diagonal cracking from the bottom left to the top right in the plan view images. The developed fracture plane has approximately the same direction as the primary orientation of the fibers, which can be qualitatively observed in the CT scans and indicates the presence of weak zones created by parallel fibers. The failure macro-crack was able to propagate through these weak zones both because the parallel fibers did not tend to span the crack and adequately restrict its growth and because crack growth parallel to the fibers may have been exacerbated by weak fiber-matrix bonding. The maximum load reached in the UC test was 823.32 kN, which was noted in the final load cycle.
The failure pattern of the DP tests (Fig. 4) exhibits half cones beneath each punch, with the half cone on the bottom of the specimen developing on the opposite side of the half cone on the top of the specimen. In addition, the two sides of the specimen in the plan view exhibit a shear differential displacement in the vertical direction, unlike what is expected for a homogenous material behavior (Fig. 1). Another discrepancy from the theoretical DP crack structure was the development of only two major radial macrocracks. These cracks can also be understood as a single diagonal crack, but they are herein technically classified as two radial cracks in keeping with traditional DP test definitions. A comparison of these diagonal crack structures with the fiber layout within the specimens provides insight into the intimate relationship between DP test performance and fiber orientation. From Fig. 4, it appears that the orientation of both the diagonal crack and the fibers in the plan view images of the DP-EXP test is largely vertical, with any horizontal component of the fiber orientations generally also running parallel to the primary failure crack. This provides further evidence that, given the similar fiber orientation layout seen in the two experiments, the DP reacted similar to the UC specimen, indicating some diagonal shear failure in a direction parallel to the prevailing fiber orientation. The maximum load reached in the DP test was 93.03 kN, which was noted in the first load cycle, and the peak load reached in the third load cycle (prior to what was taken as effective failure) was 77.78 kN, which indicates that the preceding load cycles substantially damaged and reduced the strength of the specimens.

Input
A customized algorithm (depicted in Fig. 6) was created in order to generate the distribution and orientation of fibers within a fiber-reinforced concrete specimen, based on given arrangements retrieved from the X-ray CT scanning process. Furthermore, the generated fibers that intersect the element's boundary are trimmed or removed depending on the concrete specimen production method (cored or cast, respectively). This code can accommodate any basic geometric shape for the specimens, and in this case a cylindrical specimen is assumed as per the laboratory tests. The necessary input for the program is separated in four basic parts. The first part defines the specimen dimensions. These are the diameter and the height of the matrix specimen, as well as the diameter and the length of the fiber implemented in the matrix. The second category of input is the fiber material density, for which a typical overall value for steel fibers of 7850 kg/m 3 was used. The third input value is the fiber dosage, which is expressed, again, in kg/m 3 . This input defines how heavily reinforced a specimen is, and for the cases presented herein a 177 kg/m 3 fiber content was applied. The final input determines the orientation of the fibers with respect to the spherical coordinate system, where orientations are characterized by angles h and U. The angle h represents the azimuthal angle in the x-y plane from the x-axis. The angle U represents the polar angle from the positive z-axis (cylinder axis) (see Figs. 5 and 7). As input, two distributions were imported that were fitted by appropriate distribution functions and then were used for sampling the angles of the fibers. To construct these distributions, measurements of the actual fiber orientations in the specimens were obtained from image processing of X-ray CT reconstruction images from scans prior to mechanical testing of the specimens (see Sect. 3.2 for further details). Then, this quantitative measurement data was fitted using a best-fit kernel distribution, which is a nonparametric smoothing curve representing an arbitrary probability density function, used when a continuous parametric distribution function (e.g. normal, log-normal, extreme value, etc.) cannot properly describe the data.

Generator of fiber geometry and arrangement
Based on the input noted above, the first step of the solution is to calculate the exact amount of fibers. The volume and weight of a single fiber and the volume of the entire specimen are calculated, and according to the given fiber dosage, the amount of fibers is calculated and rounded up to the closest integer. The probability density functions of the fiber orientation are derived by curve fitting the histograms of the CT investigation results. The reference axes that define the orientations  Fig. 5-right. Initially, the first point, A i , of each fiber is randomly generated inside the boundaries of the cylinder by generating a set of random coordinates (x i , y i , z i ) by means of a uniform rectangular distribution function. For the calculation of the end point, A j , of each fiber, more intermediate steps are required in order to capture the biased orientation of the fibers. The first step is to sample the angles of the fiber to be generated from the U and h orientation distributions derived from the CT investigation. The coordinates of the end point are then calculated using the sampled U' and h' angles and the predefined fiber length, which must be transformed into three dimensional cartesian coordinates (x j , y j , z j ). The result of this step of the code is the generation of preliminary fibers with input-defined length, and controlled orientation.
Following this process, the algorithm investigates whether the end-point coordinates lie outside the predefined material specimen boundaries. If the boundaries are exceeded, the fiber is rejected and a new fiber is generated by iterating the procedure defined above, starting with a new set of random coordinates for A i (x j , y j , z j ), until the conditions are met and the intended fiber dosage is maintained. This procedure mimics the rational boundary effects for a given cast sample.

Validation of generated fiber set and import to FE software
The resulting geometry was validated before importing to the FE model, based on existing data from fiber reinforced specimens that were scanned with an X-ray CT apparatus. Two validations were carried out, namely based on the statistical distribution of the fibers' angles (U', h') and based on fiber density. The results of these validations presented in Fig. 8. For validation of the fiber orientation, the same reference system is used as that from the experimental data extracted during image processing of the X-ray CT images. Orientation angles were measured from the vertical (z) axis and one of the horizontal (x) axes, and fibers have angles of 0°or 180°when they are oriented parallel to the reference axes z and x and an angle of 90°when they are oriented perpendicular to these axes (see also Fig. 5). To calculate the angle U' of each fiber, a respective vector was generated with the origin (x i , y i , z i ) and direction parallel to the cylinder axis (z), and then the angle between the fiber and the vector was calculated. For the calculation of h', the x and y cartesian coordinates of each fiber were transformed in polar coordinates and the angular coordinate h' was measured. For the density based validation, the generated field of fibers was transferred to a designing/CAD program, where it was sliced in 32 pieces with equal volume. The total length and, subsequently, the total volume of fibers in each piece is used to calculate the average fiber density in the specimen. In order to validate the accuracy of the resulting geometry, the angle distributions of the generated fibers are compared with the angle distributions obtained from the experimental data, i.e. the target distribution. This validation is conducted by overlaying the output histograms of the angles of the generated fibers on top of the angle distribution plots from the experimental results (Fig. 7). Additionally, bivariate scatters of the output and target orientation distributions are plotted, providing an overall understanding of the similarities between the experimental and numerical geometry (Fig. 8). Besides a visual comparison, which confirmed the validity of the algorithm, the Kolmogorov Smirnov test was conducted, which also confirms that the output orientation data fit with the experimentally obtained distributions. As regards the density of the fibers estimated from the CT testing results, only the mean density of fibers is compared with the experimental results in terms of overall material volume ratio. A very good agreement was also achieved between the X-ray CT measured and algorithmically generated mean fiber densities. As an overall comparison, the average fiber densities in the specimens were as follows, which shows a good agreement between the experiments and the simulations: • Unconfined compression (UC): Experiment: 3.19%; Model: 3.06% • Double punch (DP): Experiment: 3.09%; Model 3.07% After the resulting geometry is validated and representativeness of the fiber orientation in the test specimens is confirmed, the final step of the code is to export the entire geometry in the form of an IGES file containing all fibers as a linear spline element, which is imported to the FEM software. For FE modelling convenience, an intermediate step is taken to modify the spline element to individual lines (representing the fibers) by use of, for example, a commercial drawing/ CAD software. Fibers with nominal length less than ten millimeters were automatically selected and eliminated at this work-step before importing the fibers into the FE program, because such geometries can cause convergence issues and make the solution more computationally demanding without offering any additional model efficiency.

Model set-up and assumptions
For the FE simulations, the commercial software package ANSYS was used. The geometry of the models was based on the specimen geometries, and they are different in both fiber arrangement and matrix geometry. Besides the UHPC specimens, two support plates were simulated for each model, used to uniformly transfer the load via a frictional contact interface. The size of the simulated steel punches in the DP-FEM model have great importance, since the reaction force derived from the test is calculated based on the punch diameter, and because it strongly correlates with the generated fracture cone planes. The punches are simulated as plates with diameter and height of 17.5 mm. For the compression test, loading plates with diameter sufficiently greater than the cylinder diameter were simulated, i.e. with 90 mm diameter and 20 mm height. The boundary conditions are the same in both models. The bottom plate acts like a support, with its bottom face completely fixed, and the top plate acts like an actuator in both scenarios. A controlled displacement is applied on the top (i.e. external) face of the top plate along the axis of the cylinder, and the other two directions of the plate are constrained. The analyzed bodies are supported at each end by implementation of an interface with a friction coefficient of m = 0.2, which is assumed to realistically represent the experimental setup. The modelled specimens were tested under pulsating loading in line with the experiments described above. Loading of the simulated specimens included three load cycles. In the third cycle, the models were loaded beyond failure, until a clear descending load-displacement curve developed. The simulated testing protocol, which was found to most accurately reproduce the experimental behavior, is given below: • For the UC-FEM, the displacements for the 1st, 2nd, and 3rd analysis step were 0.65 mm, 0.75 mm, and up to collapse, respectively. • For the DP-FEM, the displacements for the 1st, 2nd, and 3rd analysis step were 0.40 mm, 0.41 mm, and up to collapse, respectively.
For the non-linear 3D FE analyses, the explicit analysis system of ANSYS software was used. The element used for the solid blocks was (SOLID 164), which is defined by eight nodes and has nine degrees of freedom (DOF) in each node (these comprise the kinematic definition of three dimensional nodal displacement, velocity and acceleration for explicit dynamic solvers). For the simulation of the steel fibers, a linear element (BEAM 188) was implemented, which is a linear 2-node beam with six degrees of freedom to each node based on Timoshenko theory, and the end nodes kinematically constrained with the solid elements, thus bearing both axial and transverse loads. Hence, this fiber simulation can function in bridging both tensile/flexural and shear damaged regions in the model. The reinforcement elements are simulated as fully bonded to the solid elements, on the assumption that good bonding without sliding and full utilization of the fiber is achieved, based on the literature findings discussed above in Sect. 2. This is technically achieved by constraining the nodes of the discretized fibers at their relative locations within the solid elements, and maintaining their compatibility with the solid element shape functions. Deboning is then effected for nodes in damaged solid elements, proportionally to the stiffness reduction (described in the next paragraph). Special consideration was given in the quality and the dimension of the FE mesh in the models. After a preliminary investigation, a hexahedron brick element with 8 nodes, and a maximum element size of 3 mm was selected for both models. An additional detailing of the mesh was necessary for the DP-FEM in the region of the punches, since they cause a greater stress and damage concentration toward the center of the cylinder. The refined mesh region has a diameter of 20 mm, i.e. slightly broader than the punch diameter, and an element size of 2 mm. The implemented material model maintains the same stress-strain relation independently of the mesh size and does not require a regularization algorithm when varying element sizes. With this approach, the accuracy in the main area of interest, where the cracking and failure occurs, is substantially increased. For the meshing of the fiber reinforcement, a maximum element dimension of 3 mm was chosen. The finite element meshes are shown in Fig. 9.
Damaged plasticity is currently commonly considered as an overall acceptable approach for concrete constitutive laws, and various plausible material models are available in commercial and research software. Main differentiations relate to their softening evolution equations and the degradation of the elastic stiffness. Without excluding the possibility to achieve a good representation of the tested phenomena with other modelling approaches, the Riedel-Hiermaier-Thoma (RHT) model [54] was applied for the simulation of the concrete in this study. This is an advanced plasticity model specifically for concrete and it is capable of capturing strain hardening and strain softening concrete behavior. The model assumes that damage accumulates due to inelastic deviatoric straining, which can result in stain softening and reduction in shear stiffness. The RHT model accounts for the mesh size influence by assuming a linear softening law and a proportional relationship between the element characteristic length (cube root of element volume) and the fracture energy. This can lead to some mesh dependency [55], therefore a mesh density sensitivity study was carried out prior to the analyses. The input constitutive parameters have been taken from the unreinforced material characterization tests, as shown in Table 2. All further RHT constitutive properties are calibrated by default software formulations based on the compressive strength of unconfined cubes and respective benchmarking tests [56]. It is noted herein, that for the characterization of the compressive strength of Cor-Tuf, cylinders with diameter of 75 mm were used [45], and, therefore, an additional factoring of the compressive strength by 1.1 was applied to account for the shape and size effect influence [57]. The steel fiber reinforcement was simulated as having a perfect plastic, bilinear behavior ending at plastic stain failure. The punches and the supports were simulated with a linear elastic steel material. The input material properties are listed in Table 3.

FE analyses results and comparison with CT scans and mechanical tests
The load displacement curves obtained from the FE analyses (noted as -FEM) and the laboratory tests (-EXP) are provided in Fig. 10. In the case of unconfined compression, it can be seen that the loaddisplacement curves UC-EXP and UC-FEM remain virtually linear until a load level very close to the failure load, without residual deformations or significant damage. The inclination of the load-displacement development remains essentially constant up to this point and it is virtually the same in both UC-FEM and UC-EXP. UC-FEM appears to have a more abrupt drop after the peak load as compared to UC-EXP and, hence, it appears to exhibit a smaller post-peak ductility, although it still regains a plastic behavior once it drops to 70% of the peak load. The peak loads and corresponding displacements for the UC-FEM and the UC-EXP are 767 kN/0.83 mm and 823kN/ Fig. 9 Overview of finite element model meshing: separated solid and beam elements of the UC-FEM model (above), and of the DP-FEM model (below) 0.92 mm respectively, which shows a difference of nearly 10% between both values. The double punch test response at the final loading level is influenced by the preceding loads more evidently. Although residual deformations after the first two loading loops in both DP-EXP and DP-FEM are minimal, the load reached in the initial load cycles cannot be reached in the ultimate (third) load cycle. The load-displacement response of DP-FEM is somewhat stiffer than DP-EXP, and the maximum loads reached for the same displacement are higher by approximately 15% for the DP-FEM curve. In particular, the maximum load and corresponding displacement are at 93 kN/0.45 mm and 108kN/0.4 mm for the DP-EXP and the DP-FEM, respectively.
Similarities of the numerical modelling with the experimental results are evident in the damage patterns within the specimens, as this is demonstrated through CT imaging in Fig. 4 and cross-sections from the FE models in Figs. 11 and 12. As seen in Fig. 11 (left), the UC-FEM model develops damage (highest plastic strains, in red) in a diagonal arrangement from top right to bottom left, moreover with a vertical extension toward the core of the specimen. This damage pattern is similar to the vertical cracking and columnar damage formations witnessed in the CT scans of UC-EXP in Fig. 4 (left). The stresses in the fibers in Fig. 11 (right), indicate a significant postdamage load bearing of the fibers within the core of the loaded element, and activation close to the damaged zones in the perimeter. Some fibers appear to yield within the specimen, specifically in the damaged zones. Also, most of the fibers in contact with the loading plates at the ends of the cylinder reach their yield limit locally, which is due to the fact that they initially absorb the majority of the compressive load in this region, as the stiffest material in the composite section (possibly a fictive phenomenon).
The DP-FEM also exhibits a similar damage pattern to the CT-scanned DP-EXP. In this case, the principal cracking coincides in the scanned and modelled specimen, in that its origin is located on the one side of the punch (localized load area), it propagates diagonally and it ends at the other side of the opposite punch, as seen in Fig. 12 (top). Subsequently, the model represents a rather sharp shear failure  corresponding to the one presented in Fig. 4 (right). Another significant similarity is that, instead of the typical three or four radial cracks of homogeneous materials, DP-FEM exhibits the development of two radial cracks (or, alternately, a nearly-continuous single diagonal crack), which is also the case is in DP-EXP, shown in Fig. 4 (right). As seen in Fig. 12 (bottom), the stresses in the fibers follow the damaged zones, which indicates that they contribute in the damaged regions of the specimen due to high tension concentrations. The absence of stresses in the fibers toward the periphery of the specimen (blue colored elements in Fig. 12 (bottom)), also add evidence that the fibers are inactive in non-stressed regions of the specimen, as reasonably expected. It is noted, that the deformations of the experiments presented in Fig. 10 were measured using externally mounted LVDTs (see Fig. 2). This measured displacement data could include the effects of some secondary phenomena, such as settling of the loading platens, which is also a possible underlying cause of any deviations in the modelled and measured deformations.

Conclusions
This study offers a series of insights in testing and modelling fiber reinforced high-performance materials, in this case an UHPC with a compressive strength of 237 MPa. The laboratory experiments presented demonstrate the value of X-ray CT for evaluating the internal structure of UHPC specimens, moreover in parallel to load and damaging progress. Through the use of X-ray CT, a precise visual observation of the crack patterns as well as fiber distribution and orientations is made possible. Furthermore, based on appropriate post processing of the scans, it has been made possible to obtain quantitative statistical information on the fiber orientations distribution. These data have served as an input for a high-fidelity modelling of the scanned specimens and realistic FE analyses with a discrete simulation of the fibers. The fibers were explicitly modelled as linear reinforcement elements within the concrete matrix, which was in turn simulated by three-dimensional solid elements. In order to replicate the actual orientation of the fibers in the test specimens, a customized probabilistic algorithm has been developed that samples the fiber locations and orientations derived from the CT scans. The three-dimensional nonlinear FE models have proven to adequately capture the structural response of the specimens under two complex testing configurations, i.e. pulsating compressive and double punch testing. These investigations allow for the following conclusions of the study: • X-ray CT imaging can provide not only qualitative but also quantitative fiber distribution and orientation metrics with high precision. This information can be conveyed to structural analysis models for purposes of quality assurance (regarding the material production and casting), and can be used for the creation of predictive models. • The same imaging technique, if applied with an appropriate testing configuration, can also provide a high-resolution visual representation of internal cracking processes and damage propagation inside the sample. This is valuable information for the understanding of the relationship between the material's failure performance and the fibers' orientation with respect to different failure planes. • The developed code can efficiently generate discrete fiber elements based on CT scans as a CAD file with good precision as regards the fibers' asymmetrical alignment and spatial distribution. This is verified from the generated geometry statistical validation against the input, but also based on the damage patterns obtained from FE models, once these fibers are imported to the model as structural elements. • Although certain deviations between the experimental and the numerical results are shown, given the complexity of the pulsating loading conditions and some simplifying modelling assumptions, the FE models using the explicit simulation of fibers are deemed to present a realistic load-displacement response of the tested specimens at pre-and post-damage conditions. This can enhance the understanding of the material's performance with regards to the fibers' alignment, and it can serve as a basis for inverse confirmation of material properties. This can, for example, be applied in the assessment of existing UHPC structures by use of core sampling, mechanical testing, CT imaging and identification by use of FE models. Moreover, the model can identify the stresses developed in the fibers as a reinforcement component, which can, in turn, serve as a basis for an optimization of the fiber alignment in advanced UHPC production techniques.
Besides the advantages of the material investigation approach demonstrated herein, limitations apply, and further research potential exists. As regards the model configuration, it has to be established whether and to which extend the implemented assumptions about the fibers' bond-slip characteristics and about solid matrix constitutive properties can influence the efficiency of the model, e.g. by extended laboratory testing and parametric numerical analyses. The concrete fracture performance in the vicinity of fibers and between individual fibers should be tested in a way that accounts for the size effect and complex fiber stress situations in materials with significant fiber density. Parametric analyses on the implemented fiber alignment statistical input can also disclose the model's sensitivities on this randomized input. At a further step, the proposed modelling technique may prove useful in understanding how nonhomogeneous fiber orientations, for example due to various casting methods, or due to material characterization in cast instead of cored specimens, can cause quality variabilities. Conversely, the identification of stresses in the simulated fibers can form a basis for design optimization with biased alignments for specific structural components and loading situations. In future studies, it will also be of interest to assess the effect that sources of experimental scatter (including inherent variations in mixing conditions, casting procedures, and loading configurations) have on the accuracy of the model.