Bridging Law Application to Fracture of Fiber Concrete Containing Oil Shale Ash

: Concrete is a widely used material in various industries, including hazardous waste management. At the same time, its production creates a signi ﬁ cant carbon footprint. Therefore, intensive research is being conducted to create more eco-friendly concrete, for example, partially replacing cement with by-products such as oil shale ash (OSA) or improving properties by adding dispersed ﬁ bers such as basalt ﬁ bers (BFs). The article consists of experimental testing of nine types of concrete and the modeling of crack propagation in bending. The basic trends of crack propagation in samples of concrete with OSA and BFs are simulated using a two-dimensional Finite Element (FE) model considering only material degradation on the opening crack surface and experimental data of three-and four-point bending tests. Crack propagation is modeled using the bridging law approach. A surrogate model for predicting the peak loading as a function of tensile strength and fracture work was created. An examination of the results of the FE model shows that the bilinear and nonlinear bridging law functions best describe the crack growth in the analyzed material. A comparison of experimental and modeled results showed that the length of the composite BF strongly a ﬀ ects the accuracy of the numerical model.


Introduction
Concrete is the second most used substance in the world after water [1], despite cement, one of its main compounds, being responsible for more than 7% of global carbon dioxide emissions [2]. Therefore, possibilities for producing environmentally friendly concrete have been actively studied in recent years. In 1968, Sweden was the first country to suggest to the United Nations a conference focused on human interaction with the environment, which resulted in The United Nations Conference on the Human Environment held in Stockholm in 1972. The Earth Summit in 1992, held in Rio de Janeiro, followed by the Kyoto Protocol in 1997, established a policy to protect our environment and support development. The latest international treaty on climate change is the Paris Agreement signed in 2016, which significantly accelerated the shift in human behavior and industry awareness of greenhouse gas production and resource depletion. Recent trends support the circular economy and circular material flow. In the case of concrete, the use of supplementary cementitious materials (SCM) as a cement replacement, and the use of secondary fuel resources for cement production, has significantly impacted cement's carbon footprint. The use of secondary raw materials/industrial side streams as supplementary materials fits well with the circular economy concept and helps solve the problems of climate change and environmental degradation. However, their widespread use is limited by several other factors, such as the need for more suitable applications in place and the amount of SCM locally available [3].
Adding oil shale ash (OSA) to concrete has clear advantages in reducing the use of raw materials during concrete production. At the same time, this approach helps to make good use of OSA, a byproduct of energy production generated in huge amounts due to oil shale combustion processes. The cement replacement level with fly ash is limited to 35% of the cement mass in the concrete [4]. A promising approach is using OSA instead of fly ash as an immobilizing matrix in hazardous waste management. Analysis of mineral composition and cementation properties have shown that OSA may be used as a barrier material to immobilize hazardous materials in its matrix. Studies of the leaching characteristics of concrete with OSA show that the fraction of readily soluble inorganic components decreases during curing time, and the mobility of potentially dangerous elements does not increase [5]. It should also be noted that concrete with OSA has a lower early age strength than concrete containing exclusively Ordinary Portland Cement.
The excellent mechanical properties of fiber-reinforced concrete are the reason for the growing interest in its research and application in hazardous waste's short-term and longterm storage [6]. The beneficial effects of fibers in improving post-cracking strength and energy absorption capacity are attributed to their ability to transfer tensile stress across cracks, known as the crack-bridging effect [7]. Basalt fibers (BFs) in concrete have been found to have better properties (high tensile strength, and chemical and heat resistance) and higher economic efficiency than other fibers [8][9][10]. The use of basalt fiber as a reinforcement for improving the mechanical properties of composite materials has been widely studied [11][12][13].
Concrete usually contains discontinuities such as cracks and voids that dominate the material's strength, deformability, and failure behavior. The initiation, propagation, and coalescence of cracks lead to the degradation of the mechanical behavior of concrete. Therefore, it is essential to understand the cracking behavior of concrete to ensure the reliable performance of structures [14]. Safety assessment of concrete structures used in hazardous waste management applications requires realistic models to reproduce the behavior of the materials.
The Finite Element (FE) method [15,16] has become a mainstream tool in various engineering fields, such as mechanical and civil engineering, biomedical engineering, and material science. FE is also a useful tool for modeling the mechanical properties of concrete. Khormani et al. [17] estimated the compressive strength of concrete using a 3D mesoscale finite element model created with the help of computed tomography (CT) images. The quality of the images obtained from the CT scanner allows for modeling the two types of aggregates used based on the brightness range of each aggregate. Carvalho et al. [7] used the finite element method to capture fibers' bridging effect in a four-point bending test. They use a discrete crack approach to modeling, requiring material parameters usually available in design codes and technical data sheets. A three-point bending test for high-performance concrete reinforced with short steel hook-end Dramix fibers under a 2D plane stress condition is simulated by Le et al. [18]. Fiber bridging and cohesion resistance are used as fundamental mechanisms for the development of fiberreinforced concrete models. The evolution of complex cracks, including tree-shaped cracks, regularly distributed intersecting cracks, and randomly distributed intersecting cracks, is simulated and studied by Jia [14] using the field-enriched finite element method.
Although the well-established and widely used FE method is the dominant numerical tool for solving fracture mechanics problems, other approaches are also described in the literature; for example, the boundary element method is relatively widely used [19,20]. Jiang et al. [21] demonstrate material analysis with interfacial cracks using the finite difference method. Simulations of Brazilian disc and compressive tests with the discrete element method are described in [22,23].
There are many numerical models for characterizing the initiation and propagation of cracks in concrete. One of the first nonlinear crack models to analyze the failure process of concrete was the cohesive crack model [24][25][26]. Recently, phase field theory has been increasingly used for simulations of crack initiation, propagation, merging, and branching [27][28][29]. Feng et al. [30] demonstrate the successful application of the phase-field regularized cohesive crack model to simulate the three-point bending test for concrete. On the other hand, Li et al. [31] offer to combine the phase-field regularized cohesive zone model and random aggregate models for simulating complicated mesoscale damage and fractures. Fiber-bridging constitutive law applications have been studied in [32,33].
Considering the rapidly growing availability of computing resources, there is a tendency to create complex numerical models. Our proposed approach is the oppositeto use the simplest possible numerical models to estimate the typical behavior of fiber concrete and its properties in bending tests.
In this paper, basalt fiber concrete mixes, with and without partial cement replacement by oil shale ash up to 30%, were prepared in the laboratory and characterized. A novel approach to modeling basalt fiber concrete supplemented with OSA using the discrete crack method based on the bridging law approach is presented. The bridging law curve is a material property that characterizes the macroscopic fracture behavior. The Mode I crack behavior of concrete is investigated by three-point bending (3PBT) and fourpoint bending tests (4PBT). During the tests, load-crack mouth opening displacement (P-CMOD), and load-midspan deflection (P-δ) curves are obtained. Theoretically, 4PBTs allow obtaining a pure bending condition in the beam segment between the loading rollers. Consequently, tensile cracks propagating upwards are expected to occur during the test. The P-CMOD and P-δ curves allow us to find part of the input data required for the bridging law. Four functions were tested for the materials analyzed in the study: linear [34], bilinear [35,36], non-linear [37], and exponential [38]. The most appropriate modeling results were obtained using bilinear and non-linear functions. The FE model combined with the bridging law is developed to achieve numerical simulations of 3PBT and 4PBT and obtain P-CMOD and P-δ curves, respectively. For parametric studies, it is useful to use surrogate modeling options, such as the response surface method. A data set is received for a surrogate model that predicts peak load values as a function of fracture work and tensile strength by performing a limited number of FE simulations. The flowchart of the overall study can be seen in Figure 1. The combination of low cement concrete with composite basalt fibers has not yet been investigated for behavior under bending conditions. As the hazardous waste management field is an extremely risky environment, laboratory testing does not seem to be enough. Therefore, numerical modeling was introduced to capture all potential discrepancies during experimental testing. The innovative aspect of this article is the correlation between experimental and numerical modeling of fiber concrete with OSA.
This paper is structured in the following way: The materials, specimens, and experimental tests performed are presented in Section 2. The applied numerical methods, i.e., the finite element method with the bridging law approach and surrogate modeling, are described in Section 3. All numerical examples and simulations with the proposed methodology are contained in Section 4, and the paper ends with conclusions.

Mix Design with Composite Basalt Fibers for Three-Point Bending Test Experiments
Three different mixes with and without various doses of composite basalt fibers (BFs) were designed for the experimental testing. ReforceTech supplied the basalt fibers, Mini-Bars, used in this study. Each single MiniBar composite BF is a bundle of basalt microfilaments (forming thread) impregnated by a polymeric matrix. Composite fiber's mechanical behavior in concrete was investigated by Lukasenoks A. et al. [39] and Macanovskis A. et al. [40]. Fibers show competitive properties compared with metal fibers, and interest in them is rapidly growing [41]. The length of fibers was 43 mm, double the length of the maximum aggregate size used in the mixes to ensure a bridging effect. The fibers have a diameter of 0.70 mm, density of 2000 ± 100 kg × m −3 , min value of modulus of elasticity 42 GPa, and min tensile strength of 1000 MPa [42], see Figure 2. Reference concrete mix A0 was prepared without basalt fibers. Additional mixes, A1 and A2, with a basalt fiber dose of 6.7 kg and 10.8 kg per 1 m 3 of concrete. These doses correspond to volume 0.00321 m 3 (0.32 vol%), equivalent to 25 kg of steel fibers, respectively, and volume 0.00513 m 3 (0.51 vol%), equal to 40 kg of steel fibers. Mixes were prepared with natural sand fractions 0-8 mm and crushed coarse aggregates in fractions 8-22.4 mm. The cement used in this study was blended cement, CEM II/B-S 52.5 N, with an apparent particle density of 3080 kg × m −3 . Pigment BAYFERROX 920 was used to distinguish better basalt fibers from cement matrix and, also due to additional research, focused on recycled concrete aggregate testing that is not part of this study. Water and superplasticizer Dynamon SX-23 from Mapei were used for consistency adjustments. For all cases, the W/C ratio was 0.4. The designs of analyzed cases are summarized in Table 1. Prepared concrete was tested in a fresh state for density, air content, and consistency documented by a slump. The air content measurement is essential for mixes with fibers because there is a risk of bland air by fibers. Therefore, it is an important property to monitor. The compressive strength of prepared concrete was tested on cubes with edges of 150 mm and E-modulus on cylinders with a diameter of 150 mm and height of 300 mm. Table  2 summarizes the properties of fresh and hardened (at 28 days) concrete. It can be seen from Table 2 that the lower amount of composite BF (6.7 kg per 1 m 3 ) does not significantly affect the compressive strength. In contrast, when the BF content is increased to 10.8 kg per 1 m 3 , the compressive strength decreases by 4.1 MPa (6.4%). This observation is consistent with the results available in the literature [43,44]. A test for residual flexural tensile strength is essential for fiber concrete, and, therefore, prisms with dimensions 150 × 150 × 550 mm (height × width × length) were prepared. There were 3 prisms tested for mix A0 without fibers and 6 prisms for each mix with BFs, Batches A1 and A2. All samples were prepared with a notch and tested according to NS-EN 14651:2005 + A1:2007 Test method for metallic fiber concrete-measuring the flexural tensile strength (limit of proportionality (LOP), residual) [45]. The notch was specific for the residual flexural tensile strength test, within which the maximum load value can be obtained. The notch-to-depth ratio was 0.167. The only deviation from the standard procedure was the number of tested samples, which was reduced due to technical and time constraints.
The load at the limit of proportionality was 17.4 kN for mix A0, 15.8 kN for mix A1, and 14.1 kN for mix A2. The concrete mix A1 with 6.7 kg of BFs, performed similarly to mix A2, which contains more 4.1 kg of BFs. Furthermore, the residual flexural tensile strength of mixes with BFs was comparable with the maximum reach tensile strength that was better for A1 (see Figure 3). This might be due to better fiber distribution in a sufficient amount of matrix. Higher doses of fiber can result in uneven distribution and a lack of matrix in which the fibers can be bonded.

Fabrication of Basalt Fiber Concrete with OSA for Four-Point Bending Test Experiments
Prisms from basalt fiber concrete were fabricated for mechanical experiments. Concrete ingredients were Schwenk Eco cement CEM II A-LL 42.5 N, OSA (from Eesti Thermal Power Plant, Narva, Estonia), water, dolomite gravels 4-8 mm, sand 0.3-2.5 mm, sand 0-1.0 mm, dolomite powder, plasticizer (Sika D190), and (if it was necessary according to recipe need) basalt fibers were mixed in the mixer. BFs used in the study were composite basalt fibers supplied by DBF Turbobuild INTEGRAL with dimensions 24 mm long and approx. 1 mm in diameter. According to supplier data [46], basalt rovings are coated with an alkaline-resistant material and chopped during macro fibers fabrication.
To provide conventionally made concrete, all recipes were mixed in a concrete mixer JS200L. First, dry material, cement, fine aggregates, and coarse aggregates were mixed. After mixing for 3-4 min, water (80% of the total), and plasticizer were added. While blending for another 5-6 min, water was added gradually, and concrete consistency was monitored. The selected type of mixing prevented lumps formation. After that, the control slump test was performed, and the concrete was mixed for another 2-3 min.
The procedure for a concrete prism with BF fabrication is similar; only when the consistency was ready were the BFs added to the concrete in small portions so that they did not form agglomeration and were distributed uniformly over the entire volume. After that, we poured fresh concrete into prism molds with dimensions of 100 × 100 × 400 mm (height × width × length). After 2-3 days, the prisms were removed from the molds and placed in a bath with a constant temperature of water 20 ± 2 °C. After 28 days, samples were tested in 4PBT.
All samples were marked for orientation in the bending testing machine and tested on the Automax System Controls test machine for 4PBT, according to the European standard EN 12390-5:2019 [47]. Unlike the three-point bending test, prism samples without an initial notch were tested in this case. A frame with two sensors (WA 10 mm) was placed on top of the prisms, and sensors (up to 10 mm) were placed perpendicular to the prism, so that the reading of sample deflection could be recorded. Values of the applied force and the sample deflection were obtained using the testing machine.
Depending on the dose of added BFs, the created samples could be divided into two batches-Batch B with 15.7 kg of BFs per 1 m 3 of concrete and Batch C with 23.7 kg per 1 m 3 of concrete. Mixtures with 0%, 10%, 15%, and 30% OSA were used for further analysis from the samples of Batch B. The tested compressive strength of the corresponding sample cubes with 100 mm edges was 11.7 MPa, 13.2 MPa, 11.5 MPa, and 10.5 MPa. The corresponding tensile strengths ranged from 1.13 to 1.27 MPa. Samples with 10%, 20%, and 30% OSA were used in modeling from Batch C. The compressive strength was 7.1 MPa, 7.2 MPa, and 6.6 MPa, respectively. The tensile strength ranged from 0.93 to 1.3 MPa. These results do not coincide with common views about increasing the concrete strength if you add fly ash or OSA in small amounts. It is necessary to point out that we are investigating quite high OSA concentrations.
The results show that the compressive strength of samples with OSA also decreases when the dosage of BFs is increased. By increasing the number of fibers by 1.51 times, the compressive strength decreases from 1.86 to 1.59 times, depending on the amount of added OSA. The reason for this is the higher fiber volume fraction in the concrete mix, which is more challenging to mix uniformly. The workability of mixes in Batch B tested by slump test was between 20 and 30 mm. The addition of 10% OSA to the concrete mix increased the compressive strength by 12.5%, but further addition of OSA had the opposite effect. For Batch C samples, adding 10% or 20% OSA to the concrete mix practically does not affect the compressive strength, but further increasing the amount of OSA decreases the compressive strength value. The E-modulus of the studied concrete with BFs and OSA ranged from 16 to 30 GPa. The typical modulus of elasticity for basalt fibers is 87 GPa [44]. The structural steel from which the rollers of the bending test machine were made has an E-modulus of 210 GPa. From the above, it can be assumed that there were no noticeable deformations in the loading and supporting rollers of the test machine, and they can be ignored when modeling the bending test.
Of all the analyzed BF concrete mixes with added OSA, the highest compressive strength value was obtained for the mixture composition with 15.7 kg of BFs per 1 m 3 of concrete and 10% OSA. However, considering the insufficiency of the statistical amount of the analyzed experiments, it is not yet possible to determine the optimal composition of BFs and OSA in concrete at this research stage.

Finite Element Model
Using the Finite Element (FE) method, 3PBT and 4PBT were modeled to evaluate the typical crack propagation in BF concrete samples with and without OSA in its mixture. The bridged crack propagation was modeled by adapting the in-house software described in [48]. A detailed description of the software can be found in Appendix A. Analyses were made of both the crack propagation path and the crack opening, which is a non-trivial function of the crack path as a function of loading.
The propagation of a crack in the prism beam was numerically modeled using the principles of force balance and energy minimization. Several assumptions were introduced in the modeling, which helps to evaluate the trends of changes in the mechanical properties depending on the parameters, even if they do not correspond to the reality for each specific sample. A small displacement/angle approximation was used in the FE model, as experiments show that the angle formed between the crack surfaces is small. The results were obtained without deforming the computational mesh. Material degradation in volume was not modeled. Significant material degradation was assumed to occur only on the surface. The material was considered to be isotropic with a given modulus of elasticity and Poisson's ratio. The presence of a notch in the sample was considered in modeling the 3PBT. In the 4PBT, it was assumed that the crack will start forming in the middle of the prism in the longitudinal direction (x-axis); see Figure 4a. In the 3PBT and 4PBT simulations, the crack formed two symmetrical surfaces with respect to the y-axis (see Figure 4b). The crack propagated in a straight line upwards; therefore, only half of the prism sample was modeled. The results described in [49] show that although the crack tends to propagate along weak points and thus become tortuous due to random inhomogeneity, such an assumption is acceptable when we are interested in the "average" typical crack path. The linear-elastic 2D problem, chosen for its simplicity, has inherited well-known limitations, which can be divided into three groups. 1. Two-dimensional assumption of flat, straight crack lips in sectional plane YZ at the geometrical center of the specimen rather disagrees with the experimental observations. Due to the stochastic nature of concrete's microstructure, the overall crack propagation mechanisms have evident 3D traits. 2. The fiber pullout and microcracking contribute to the volumetric degradation of the material near the crack path, which is unaccounted for by the linear elastic model. 3. Finally, all deformations and displacements are considered small. Any finite solid body rotation during the load is neglected, which may play some role in the final phase of the cracking. Some of these inaccuracies (1-3) can be addressed in further work. The authors have already considered the 3D probabilistic approach instead of the 2D determinism used so far. However, this particular investigation was driven by practical considerations, and the authors were cautious about the possible rapid increase in the model's complexity in an expansion of the principal simplicity. Despite potential improvements in accuracy, more complex models may be even less favorable in the engineering practice-notable examples are difficulties and expenses related to the identification of additional material properties. Overall implementation availability in an FE analysis code also must be counted on. The higher probability of human errors with the application of conceptually more complex, sometimes even obscure, models includes all of the abovementioned, and it has to be considered in a further model's refinement.
The crack opening condition is imposed on the edge AB ( Figure 4). Edge CD is free to move in the XY plane. No bridging law is applied to the upper crack path point (see point B in Figure 4b). The modeled bridging law transfers tensile load in a single, x-axis direction. Thus, once the crack path is opened at the particular point (node), only normal traction (positive normal direction) forces can be transferred over the crack faces. One can consider the integral force balance for the whole region (½ of the specimen). Since no other loads act at the x-direction and only traction in the x-direction is allowed along the bridged crack part, some material volume must be under compression (σxx < 0) on the crack path near the opposite edge (see Figure 4). While the crack grows, the compressive region's length shrinks, and the compressible stress component's absolute value increases. It is known [48] that the material may fail under compression at certain bridged crack sizes in polymeric composites. Still, preliminary analysis shows that this does not happen in the concrete of the type analyzed here due to large differences in tensile and compressive strength. In the model, the need for the force balance mentioned above is considered by fixing a single node (point B) with no crack opening at the final stages when the rest of the crack path is already abridged.
The external load, defined by the displacement in time, is applied to the contact surface of the loading roller. Loading and supporting rollers are modeled as contact surfaces (see A'D' in Figure 5), to which a specific algorithm for energy minimization is applied. When modeling bending tests, it is essential to adequately describe the interaction of the supporting and loading rollers with the specimen. The contact points between the sample and the rollers vary according to horizontal and vertical displacements. Mujika et al.'s [50] study considered the effect of changes in the contact point between supporting and loading rollers in three-point and four-point bending tests and the impact of possible horizontal reactions when small bending angles are assumed. They conclude that the effects of contact point displacements and horizontal forces are more critical in four-point bending than in three-point bending. Energy, introducing the generalized coordinate q, can be formulated as follows: where the K matrix is symmetric, and qi can be either a displacement or an angle. The corresponding matrix K constants are found using the FE model, which solves the following six sub-models in each iteration: 1. The constant displacement du s = 0 and the rotation angle θ s = 0 is defined for the support (lower) roller. The roller of the test machine (upper) is loaded with the displacement du t , but the angle θ t = 0 (see Figure 4). 2. The constant displacement du s = 0 is defined for the support roller; it is loaded with the rotation angle dθ s . For the roller of the test machine, the displacement is constant du t = 0, and the angle is zero. 3. The support roller has the constant displacement du s = 0 and the rotation angle θ s = 0, while the roller of the test machine is loaded with the rotation angle dθ t , and the displacement is constant du t = 0. 4. The support roller has the constant displacement du s = 0, but it is loaded with the rotation angle dθ s , while the roller of the test machine is loaded with the displacement du t , but θ t = 0. 5. The support roller has a constant displacement du s = 0, and the rotation angle is zero.
The roller of the test machine is loaded with both the displacement du t and the angle of rotation dθ t . 6. Both rollers are loaded with the angle of rotation, but the displacements of both rollers are constant.
After finding the constants, energy minimization is performed. Knowing the energy (Equation (1)), we can apply the second theorem of Castigliano to determine the forces by Equation (2) and the moment by Equation (3).
where u is the displacement, and θ is the angle, see Figure 4b.

Analyzed Bridging Law Functions
The bridging law is the main input required by the FE model to characterize both the matrix and the fibers. Parameters such as the effects of moisture ingress or environmental aging, as well as fiber orientation, volume fraction, or surface treatment, are not defined in the FE model, but by the chosen bridging law. Comparing the modeling results with the experiments, bilinear and nonlinear bridging law functions were used for further modeling. Both functions require the following parameters to characterize them: the tensile strength ft and the fracture energy GF.
Based on Herrmann and Roux [51], the fracture energy GF is defined as the area average of the fracture work calculated by Equation (4).
where W represents the work for generating the new surface, and Alig is the ligament area formed at the front of the initial crack. The fracture work absorbed by the material throughout the fracture process in the bending test can be calculated by Equation (5).
where P represents the external loading, du represents the change in the CMOD (or δ), and a and b are the lower and upper limits of the CMOD (or δ) [52].

Bilinear Bridging Law Function
The expressions for the bilinear bridging law function described by Guinea et al. [35] are the following: where σ is stress, w1 = 1.2 (GF/ft) is the horizontal intercept of the steeper segment, and σk, wk are the coordinates of the kink point, given by Equations (7) and (8), respectively.
The critical opening wc is obtained from the quadratic equation: where w is the abscissa of the centroid of the bridging law curve calculated by Equation (10).
α is a dimensionless parameter, which, for a bilinear function, is 0.987 according to Petersson [53] or 1.1-1.3 according to Reinhardt et al. [54]. In the analyzed cases, the best agreement with the experiment was given by α = 1.1.

Nonlinear Bridging Law Function
The nonlinear bridging law function developed by Cornelissen et al. [37] is proposed in the following form: where c1, c2, and wc are material constants. For analyzed concrete, the three parameters are considered as c1 = 2.5, c2 = 6.93 and the critical opening is obtained by Equation (12).

Surrogate Modeling
The use of surrogate models helps to effectively understand and evaluate the interdependence of model factors (inputs) and responses (output) without full-scale numerical simulations or experiments. In the general case, the goal is to obtain a response function F (as in Equation (13)) that depends on the values of the factors xi.
where x1, x2, … xn are coordinates from the factor space. The response function F in the factor space represents the response surface. This research examined the first, second, and third-order polynomial approximations. The coefficient of determination R 2 and Pearson's χ 2 test were used to characterize the fit of the surrogate model. In addition, a leave-one-out cross-validation approach was used for surrogate model verification. Leave-one-out cross-validation is a cross-validation method that involves using one observation as validation data, and the remaining data are used to train the model. For example, point i = 1 is selected for validation in the first iteration for a data set with N points, and the remaining (N-1) points are used for model building. For the second iteration, point i = 2 is selected for validation, and the remaining ones are for building the model. The process is repeated N times. A smaller cross-validation value means a higher accuracy of the surrogate model. The cross-validation percentage error is defined in Equation (14) to avoid overestimating the prediction error.
where is the approximation, ignoring point i (i = 1, 2, ..., N). A surrogate model is assumed to have acceptable accuracy if the cross-validation error is less than 10% for numerical experiments and less than 5% for physical experiments.

Numerical Examples
To validate the suitability of the FE modeling approach discussed above for studying the mechanical properties of fiber concrete, modeling was carried out for 3PBT and 4PBT using the same sample batches as in the experiments. Considering that prisms with a single notch are used in the 3PBT experiments, the nodes of the corresponding elements on both sides of the crack patch are also marked as unconnected in the FE model before calculations. The existence of a notch fits well with the model's assumption that the crack will originate in the middle of the sample in a horizontal direction and propagate vertically upwards. Three notched concrete prism examples without and with BFs under 3PBT and examples of concrete beams with two doses of BFs and three doses of OSA under 4PBT were simulated. A surrogate model describing the maximum load Pmax as a function of the tensile strength ft and the fracture work W was created using the 4PBT modeling results for a concrete specimen with BFs and OSA.
FE models with corresponding dimensions are shown in Figure 5 for both bending tests.
To evaluate the influence of the size of the cell elements on the obtained result, three types of meshes were examined for the concrete sample from Batch C. The calculated maximum force values depending on the maximum size of mesh elements are summarized in Table 3. Considering the obtained results and the time required for computing, a mesh with 2456 elements was chosen as the optimal case for further calculations. An example of a computational domain mesh can be seen in Figure 6. In addition, the refinement of mesh elements was applied to the crack path AB (see Figure 5) and the contact surfaces of the rollers with the sample (see A'D' in Figure 6). First-order triangular finite elements are used for simulations. Linear triangular elements have been widely used in engineering problems because of their adaptability to complex models, simple calculations, and acceptable accuracy [55].  The display of the bending test result-P-CMOD or P-δ curve can be divided into three parts, see Figure 7. The first ascending part of the curve corresponds to the linear phase. After the linear part, the matrix failure phase begins, in which cracks start to propagate in the concrete matrix. The load reaches its maximum value. In the third phase, continuing to increase the external load, the fibers begin to work-they are pulled out and/or broken. An example of the stress distribution in the Batch C sample obtained by the FE model is shown in Figure 8. Figure 8a shows that the solution in the linear phase corresponds to the Euler-Bernoulli beam theory. The upper left corner of the sample shows compression, while the lower left corner shows tension. Places of loading and supporting rollers are also well visible. While in the matrix failure phase, when the crack starts to open, the stress field has changed significantly, see Figure 8b. The tension-dominated area of the crack path has increased considerably.

Three-Point Bending Tests
The P-CMOD curve obtained in FE modeling and measured in experiments for samples in the three-point bending test can be seen in Figure 9. The reference samples are without basalt fibers (Batch A0); see Figure 9a. In addition, samples with 6.7 kg per 1 m 3 (see Figure 9b The maximum load of the concrete sample with 6.7 kg per 1 m 3 of BFs (Batch A1) is 15.8 kN when CMOD = 0.036 mm. The comparison of the experimental and modeling results of Batch A1 shows that the bilinear bridging law function better describes the trends of the linear part of the P-CMOD curve and cracking phase in the concrete matrix but incorrectly predicts the fiber failure phase. The nonlinear function overestimates the peak load (Pmax = 17.76 kN at 0.054 mm). However, it best describes the fiber failure phase.
Adding 10.8 kg of BFs per 1 m 3 to the mixture (Batch A2), tested on six samples, reduces the peak load value to 14.1 kN. CMOD value is 0.033 mm. It can be seen in Figure  9c that the FE model predicts a significantly higher load peak value (Pmax = 17.32 kN at CMOD = 0.074 mm) and a slower sample rupture than was obtained in the experiment. The addition of longer composite BFs causes considerable damage to the specimen volume in the bending test. Therefore, the assumption used in the FE modeling about damage formation only around the fiber surface is not true in this case. For this reason, a faster decrease in the P-CMOD curve is observed in the experiment after the maximum loading than predicted by the FE model. Modeling such cases (with longer fiber pull-outs) should also consider micromechanics models.
The results of the FE model of the 3PBT show that for the bending test modeling of BF concrete samples, a more accurate P-CMOD curve is obtained using the nonlinear bridging law function, and for plain concrete (0.0% BF)-a bilinear function.

Four-Point Bending Tests
The comparison of the experimental and FE model results of the 4PBT shows that the bilinear bridging law function best fits the experimental results of Batch B (the dose of BFs is 15.7 kg per 1 m 3 ). It follows from Table 4 that the maximum load obtained with the FE model is close to the experimentally measured values. The standard deviation stdev of four experiments is shown in brackets. The midspan deflection δPmax, at which these maximum values are reached, is 0.05 mm in the experiments, regardless of the amount of added OSA, but in the FE model, the deflection is larger, i.e., 0.09-0.1 mm. However, comparing the P-δ curve obtained in the experiments with the results of the FE model, it can be concluded that the model predicts the general trends of crack propagation. In the case of Batch C samples (the dose of BFs is 23.7 kg per 1 m 3 ), the experimental P-δ curve is most accurately represented by the FE model with a nonlinear bridging law function. The obtained results are summarized in Table 5. It can be seen that the FE model can predict the maximum loads and the corresponding midspan deflection values of the samples of this group according to the experimentally measured values under the 4PBT. Summarizing the obtained modeling results of 4PBTs for samples with BFs and OSA, the bilinear function is the most suitable bridging law for samples with a smaller number of added fibers. On the other hand, if the BF dose is 23.7 kg per 1 m 3 , a more accurate P-δ curve can be obtained with a nonlinear bridging law function. At higher BF doses, the model fits well with the experimental data for the peak loading Pmax and the corresponding midspan deflection δPmax value. Although preliminary results show that the nonlinear function gives a more accurate bending test result in mixtures with higher fiber content, more research is needed to confirm this conclusion.
The results show that adding OSA to the mixture of Batch B samples does not significantly affect the maximum load values. On the other hand, the load peak value for Batch C samples increases with the addition of OSA; i.e., with 10% OSA, the maximum is 7.44 kN, while with 30% OSA, it increases to 9.64 kN.

Surrogate Model to Predict Maximum Load in 4PBT
The 13-point orthogonal design of experiments [56] is chosen to create a surrogate model for predicting maximum loads depending on two factors (tensile strength ft and fracture work W) for the 4PBT. In the design of experiments, orthogonality is essential because it leads to the best models and reduces the variance of the regression coefficients. The building set uses BF concrete supplemented with OSA, as in Batch C. The design of the experiments is shown in Figure 10. Based on the experiments, the following factor space was chosen: ft = [0.98; 1.31] (MPa) and W = [6.5; 12.7] (Nm). The validation set consists of five experimental data points of Batch C. The cross-validation error and the coefficient of determination of the analyzed three polynomial approximations are summarized in Table 6. The coefficient of determination R 2 for all cases shows that the surrogate model agrees with the FE model. On the other hand, from the calculated cross-validation error, it can be concluded that the first-order polynomial approximation provides the worst fitting to the FE model results, and the third-order polynomial approximation provides the best. The cross-validation error of the surrogate model using a third-order polynomial for approximation is less than 0.5%. In all cases, Pearson's χ 2 criterion is met.
The response surface obtained by the third-order polynomial approximation (Equation (15)) is shown in Figure 11. The experimental points of the 4PBT are shown with green dots. It can be seen from Figure 11 that in certain areas of the design space, the Pmax values predicted by the surrogate model are lower than the experimentally measured ones. When calculated, the coefficient of determination is 0.764. To increase the accuracy of the surrogate model, the methodology for finding the tensile strength ft should be improved.

Conclusions
The study focused on the correlation between the experimental and numerical simulation of a new type of concrete containing various doses of oil shale ash and two types of composite basalt fibers. In total, nine concrete mixes were designed and tested in a fresh and hardened state. Batch A contained BFs with a length of 43 mm and was characterized by E-modulus between 28.7 and 29.8 GPa and residual flexural tensile strength test with loads at the limit of proportionality between 14.1 and 17.4 kN. Batches B and C were designed with various doses of BFs with a length of 25 mm and a replacement of 10, 15, 20, and 30% of cement by oil shale ash. The maximum load of four-point bending tests ranged between 7.44 and 10.45 kN.
Experimental results from three-point and four-point bending tests were compared to numerical simulations. In developing the finite element model, an approximation of small angles and deformations was adopted, corresponding to the experimentally observed behavior. The modeling considered changes in the contact surfaces between the concrete prism and the rollers during the bending test.
For studying the properties of new materials (with several parameters), numerical modeling provides obvious advantages compared to experiments-preparing samples and performing mechanical tests take weeks (tests are typically performed after 28 days) or longer (after 56 or 90 days), while obtaining modeling results takes a few hours.
The most appropriate bridging law is a bilinear function for samples with or without a small dose of shorter basalt fibers (24 mm). On the other hand, if the BF dose is higher, a more accurate P-δ curve can be obtained with a nonlinear bridging law function. However, it should be noted that these are preliminary results that need to be verified using an extensive set of experimental data.
Considering that adding longer composite BFs (43 mm) causes significant damage in the sample volume during the bending test, the applied FE modeling methodology can only provide a general assessment of the bending test parameters.
A surrogate model for predicting the maximum load as a function of tensile strength and fracture work was created. From the calculated cross-validation error, it can be concluded that the third-order polynomial approximation best fits the FE modeling results. The cross-validation error was less than 0.5%.
In further research, the authors would like to focus firstly on collecting more experimental data for reliable calibration of the models and applying a probabilistic description of cracking with the involvement of 3D models with finite displacements for bridging-law identification while leaving the 2D approach still valid for structural optimization. The second priority is the description of crack-related volumetric degradation of the material.  Data Availability Statement: Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A
The software used is the home brew research Finite Element Method (FEM) code used by one of the authors for the last two decades [57]. The backbone of the code is written in plain ANSI C language for POSIX (UNIX/Linux) type OS and consists of ~100,000 lines of code. The code can be grouped into the following blocks (only solid mechanics related): (1) Data input/output and memory management; (2) Two-dimensional mesh generator and 3D mesh automatic generator by mesh extrude method; (3) Common sparse FEM matrix generator by energy minimum and Galerkin methods; (4) Solver of linear real/complex equation system, direct, and Krylov space methods; (5) Plastic flow modulus, dynamic, Prager, and Armstrong-Frederich types of plasticity; (6) Linear-elastic, thermal-elastic and plastic/triangular, tetrahedral, and beam elements; (7) The procedural script interpretation modulus gives access to all numerical modules mentioned, used for boundary conditions and extra physical feature coding.
The abovementioned scripting facility is used for bridging-law-related objects and loading sequence implementation. A typical script modulus is an ASCII file with fixed command sequence (more than 150 internal commands) available as well as real, complex, integer, character data types, arrays, and matrix operations.
During the years, the code has been intensively cross-validated with the commercial code ANSYS (registered trademark of ANSYS, Inc, https://www.ansys.com/ accessed on 15 June 2023) (by direct comparison of stress-strain fields and integrally) as well as with some well-known linear-elastic analytic solutions (like a plate with a hole, Euler beam, etc.). The convergence and stability of iterative and direct solvers are controlled in code by residual analysis. Eigenvector/eigenvalue analysis of the FEA matrix can also be used; however, it is usually not necessarily due to the dominance of the diagonal terms for the linear-elastic problem on well-generated Delaunay meshes [58].
The script code fragment below calculates nodal forces: newprg calcforce