Analyzing the Validity of Brazilian Testing Using Digital Image Correlation and Numerical Simulation Techniques

: Characterizing the mechanical behavior of rocks plays a crucial role to optimize the fracturing process in unconventional reservoirs. However, due to the intrinsic anisotropy and heterogeneity in unconventional resources, fracture process prediction remains the most significant challenge for sustainable and economic hydrocarbon production. During the deformation tracking under compression, deploying conventional methods (strain gauge, extensometer, etc.) is insufficient to measure the deformation since the physical attachment of the device is restricted to the size of the sample, monitoring limited point-wise deformation, producing difficulties in data retrieval, and a tendency to lose track in failure points, etc. Where conventional methods are limited, the application of digital image correlation (DIC) provides detailed and additional information of strain evolution and fracture patterns under loading. DIC is an image-based optical method that records an object with a camera and monitors the random contrast speckle pattern painted on the facing surface of the specimen. To overcome the existing limitations, this paper presents numerical modeling of Brazilian disc tests under quasi-static conditions to understand the full-field deformation behaviors and finally, it is validated by DIC. As the direct tensile test has limitations in sample preparation and test execution, the Brazilian testing principle is commonly used to evaluate indirectly the tensile strength of rocks. The two-dimensional numerical model was built to predict the stress distribution and full-field deformation on Brazilian disc under compression based on the assumptions of a homogenous, isotropic and linear elastic material. The uniaxial compression test was conducted using the DIC technique to determine the elastic properties of Spider Berea sandstone, which were used as inputs for the simulation model. The model was verified by the analytical solution and compared with the digital image correlation. The numerical simulation results showed that the solutions matched reasonably with the analytical solutions where the maximum deviation of stress distribution was obtained as 14.59%. The strain evolution (normal and shear strains) and displacements along the central horizontal and vertical planes were investigated in three distinguishable percentages of peak loads (20%, 40%, and 90%) to understand the deformation behaviors in rock. The simulation results demonstrated that the strain evolution contours consistently matched with DIC generated contours with a reasonable agreement. The changes in displacement along the central horizontal and vertical planes showed that numerical simulation and DIC generated experimental results were repeatable and matched closely. In terms of validation, Brazilian testing to measure the indirect tensile strength of rocks is still an issue of debate. The numerical model of fracture propagation supported by digital image correlation from this study can be used to explain the fracturing process in the homogeneous material and can be extended to non-homogeneous cases by incorporating heterogeneity, which is essential for rock mechanics field applications.


Introduction
Characterizing the mechanical behavior of rocks plays a key role in optimizing the fracturing process and enhancing production from unconventional reservoirs. However, due to complex mechanical stratigraphy and heterogeneity, sustainable economic oil and gas production is still the most challenging issue as a consistent prediction of the fracture process is critical [1][2][3]. Mitigating these problems requires the accurate determination of mechanical properties such as stress magnitudes, load-deformation nature, full-field deformation, and reasonable rock strength [4]. Therefore, an in-depth understanding of the deformation characteristics of rock is essential for sound recovery from unconventional resources.
The mechanism surrounding how a rock failure occurs is a complex process of several distinct intervals including stress accumulation, interaction, initiation of fracture, propagation, and final failure [5]. Although researchers had been trying to solve the problem of the complex rock failure process using local measurement systems (strain gauges, extensometer, LVDT, etc.), these were not accurate and sufficient enough for modeling the deformation behavior of rocks. Due to the sample preparation and installation between a rock specimen and loading platens, the measurement of rock deformation and tensile strength by direct tensile testing is a challenging task [6].
Brazilian testing is one of the most widely used indirect tensile testing methods mainly due to its simplicity of operation. In this test, a circular and thin disc specimen (concrete, rock, and rock-like material) is diametrically loaded [7] to obtain tensile strength [8,9].
Numerous experimental and numerical studies had been conducted to analyze stress distribution and strength determination during Brazilian testing [10][11][12][13][14]. However, the Brazilian testing under-estimated the strength of the rock; besides, in terms of failure mode and fracture initiation point, the validity of this test is still an issue of debate. According to the findings of some researchers [15,16], the center of the Brazilian disc is the fracture initiation point, which was substantiated by experimental and numerical investigations. However, the investigations of some other researchers [17][18][19][20] supported that the crack initiation point was near to the loading point, away from the center. This inconsistency between both theoretical and laboratory analysis is likely to be due to the assumptions, criteria, and heterogeneity of the tested rocks. However, for characterizing the deformation behaviors, it is essential to accurately determine the mechanical properties of rocks. It is particularly crucial to gain a thorough understanding of fracture propagation in both laminated and anisotropic rocks.
When researchers had been trying to solve these complex failure processes using traditional methods, optical techniques provided several additional advantages [21][22][23], such as identifying fracture generation and propagation, strain measurement, and monitoring of rock deformation under loading condition.
Digital Image Correlation (DIC) is a very powerful non-contact measurement technique that is used for capturing full-field and real-time deformation behaviors for any surface [24][25][26]. This technique introduced by Chu et al. [27] established a widely spread application in material testing [28]. Mokhtari et al. [22] applied this method to characterize the fracture initiation and propagation on naturally fractured rock. Similarly, application of DIC was extended to characterize the deformation behaviors in sandstone and carbonates rocks [29,30], laminated shale [23,31], Buda limestone [32], and Mancos shale [33].
There was significant attention by the researchers on rock failure instead of predicting the pre-failure stress distribution over time. Some researchers claimed that the initiation point of failure was decided by the ultimate loading arc angle in their numerical [34,35]. Based on the homogeneous, isotropic, and linearly elastic numerical model, more numerical studies have considered the effects of anisotropy and heterogeneity on tensile strength [34,[36][37][38]. However, fracture initiation point Brazilian testing is still an issue of debate. Although anisotropy and heterogeneity has caused the change of stress distribution pattern, the verification of the accuracy of the homogeneous and isotropic numerical model as a benchmark has played a more critical role in locating the initiation point of failure. Some single-point comparison of measurements has been made in the literature, but limited Energies 2020, 13, 1441 3 of 18 studies have explained the detailed comparison of full-field deformation at various stages of loading. Moreover, to measure the deformation and tensile strength of anisotropic rocks, the formulas were less developed due to the difficulty of experimental validation [39,40]. This required more experimental and numerical studies to accurately evaluate the geomechanical properties of rocks and evaluate factors impacting rock deformation.
In this paper, we built a numerical model using COMSOL Multiphysics to predict the stress and strain distribution and displacement evolution before the failure occurred for a homogeneous rock as a benchmark. To get the elastic properties of rock for simulation input, a uniaxial compression test was conducted in Spider Berea sandstone. Brazilian testing was conducted using DIC to monitor the deformation on the entire surface of the specimens function of time. Experimentally obtained strain evolution (tension, compression, and shear) and displacement changed along the central horizontal and vertical planes, which were investigated at three distinguishable percentages of peak loads (20%, 40%, and 90%) to understand the deformation behaviors in rock. The simulation model results at these three distinct loads were compared with 3D DIC results to investigate the validity and limitations of the linear elastic homogeneous simulation model. It was found that the model could be used to predict stress distribution, full-field deformation, and fracture behaviors for different types of rocks by changing the relevant elastic parameters.

Sample Preparation
To evaluate the deformation behavior of the rocks in this study, we performed indirect tensile and uniaxial compression testing. The homogeneous and isotropic Spider Berea sandstone was collected from Kocurek Industries [41]. Then, the core samples were cut into Brazilian discs using a diamond blade to meet ISRM [42] and ASTM [7] standards dimensions (2" Diameter × 0.6" Thickness) from long cores (2" Diameter × 10" Thickness).
To determine the rock elastic properties (Young's modulus and Poisson's ratio) for simulation inputs, the uniaxial compression test was conducted in the sandstone sample (about 2" Diameter × 2.2" Thickness) following the ASTM D 2938 [7]. We cut the sandstone cores and polished it using a wheel (7" Diameter) made with nickel-coated diamond grit and synthetic resin bond to get the finest flat ends ( Figure 1) and meet ISRM [42], ASTM D 3967 [43], and ASTM D 2938 [7]. The sample properties are listed in Table 1 [41]. To remove all pore fluid, the sample was kept in a vacuum oven at 115 F for 24 h.    Spraying the paint on the front face of the sample was the final step in specimen preparation, which resulted in the black spot on the white base speckled pattern, which was prepared using Rust-Oleum flat paint ( Figure 2). The image correlation software monitored the displacement change of the applied speckled pattern and the software used the measured displacement of the speckles to calculate the full-field strain map of the sample. Generally, the local stress-strain distribution might be impacted by coating [44] while spraying on the surface of the sample disc but the paint-coated layers were maintained as thin as possible (<100 µm) to minimize the effect [45]. The thin layer ensures the paint moves with the base rock material and won't affect the displacement behavior. Therefore, the effect of the paint was considered negligible on the displacement properties of the sample. Nath and Mokhtari [23] described the procedure for creating the speckle pattern in the rock sample face.

Deformation Measurement Using a Digital Image Correlation Method (DIC)
The DIC is an optical method used to determine the displacement of an object [25] by examining the pixel intensity of speckled images before and after deformation. The DIC system consists of an image acquisition system and correlation software.
At the beginning of the analysis, a reference image was taken where a region of interest was defined in the specimen surface by dividing it into the numbers of subsets. The displacements were calculated at each and every point of the created subsets from which the full-field deformation was obtained. The deformation tracking process was reached using the search algorithm corresponding to the pre-defined similarity criterion for reference and deformed images. When the deformation of both reference and target subset centers was obtained, the maximum similarity was attained in the deformed image. For details of the DIC technique and working principle, interested readers can explore the relevant literature of different researchers [21,23,25].

Deformation Measurement Using a Digital Image Correlation Method (DIC)
The DIC is an optical method used to determine the displacement of an object [25] by examining the pixel intensity of speckled images before and after deformation. The DIC system consists of an image acquisition system and correlation software.
At the beginning of the analysis, a reference image was taken where a region of interest was defined in the specimen surface by dividing it into the numbers of subsets. The displacements were calculated at each and every point of the created subsets from which the full-field deformation was obtained. The deformation tracking process was reached using the search algorithm corresponding to the pre-defined similarity criterion for reference and deformed images. When the deformation of both reference and target subset centers was obtained, the maximum similarity was attained in the deformed image. For details of the DIC technique and working principle, interested readers can explore the relevant literature of different researchers [21,23,25].

Experimental Procedure
The experimental procedure for Brazilian testing with DIC is shown in Figure 3. The Spider Berea sandstone was placed in a Brazilian test frame [42] showing the front side of the painted sample. The Instron load frame was used to apply the diametrically opposed compressive loads subsequently while in a synchronized way, the camera was tracking the image of the speckled face of the sample. The procedure for the camera was sets and calibration as described in the ARAMIS Professional software manual [46]. Two cameras were connected to an ARAMIS computer for storing the images and operating the DIC software. The stored images were used to analyze the deformation behavior using ARAMIS 3D DIC software.
Energies 2020, 13, x FOR PEER REVIEW 5 of 20 The experimental procedure for Brazilian testing with DIC is shown in Figure 3. The Spider Berea sandstone was placed in a Brazilian test frame [42] showing the front side of the painted sample. The Instron load frame was used to apply the diametrically opposed compressive loads subsequently while in a synchronized way, the camera was tracking the image of the speckled face of the sample. The procedure for the camera was sets and calibration as described in the ARAMIS Professional software manual [46]. Two cameras were connected to an ARAMIS computer for storing the images and operating the DIC software. The stored images were used to analyze the deformation behavior using ARAMIS 3D DIC software. Firstly, at no-load condition, the reference images were taken for calibration purposes. After checking the calibration parameters, the cameras were fixed to capture two images per second during the entire test keeping at a constant displacement rate of 0.05 mm/min. Two blue LED lights were positioned in front of the Instron load frame. The lights were kept at an equal height with the sample to make sure the specimen speckled pattern was clearly recognized by the DIC camera. According to the ASTM D 3967 [43], the specimen fails in about 5 min, which fell within the suggested time frame of 1-10 min.
During the Brazilian experiment, two digital cameras (12 MP) captured an image size of 4096 × 3000 pixels that were used to monitor the field of view, which was 76 mm × 56 mm. The typical noise level in DIC strain measurement is 0.01% and the displacement error is 0.01 pixels [46]. The DIC analysis was conducted using a subset size of 30 pixels and a point distance of 15 pixels. We got a displacement accuracy of 0.2 μm in our experiment.

Brazilian Testing
The Brazilian tensile strength for a cylindrical sample is measured by the relationship [7], Firstly, at no-load condition, the reference images were taken for calibration purposes. After checking the calibration parameters, the cameras were fixed to capture two images per second during the entire test keeping at a constant displacement rate of 0.05 mm/min. Two blue LED lights were positioned in front of the Instron load frame. The lights were kept at an equal height with the sample to make sure the specimen speckled pattern was clearly recognized by the DIC camera. According to the ASTM D 3967 [43], the specimen fails in about 5 min, which fell within the suggested time frame of 1-10 min.
During the Brazilian experiment, two digital cameras (12 MP) captured an image size of 4096 × 3000 pixels that were used to monitor the field of view, which was 76 mm × 56 mm. The typical noise level in DIC strain measurement is 0.01% and the displacement error is 0.01 pixels [46]. The DIC analysis was conducted using a subset size of 30 pixels and a point distance of 15 pixels. We got a displacement accuracy of 0.2 µm in our experiment.

Brazilian Testing
The Brazilian tensile strength for a cylindrical sample is measured by the relationship [7], Energies 2020, 13, 1441 where BTS is the tensile strength of the Brazilian disc, the load P max is the diametrically opposed maximum load, D is the disc diameter, and l is thickness. Equation (1) is generally used for homogenous media and the formula is generally used for strength comparison of various rocks. For any point (x, y) in the sample in a Cartesian plane (Figure 4b), the analytical stresses solutions in the x-direction (σ x ) and y-direction (σ y ), and shear stress, (τ xy ) can be expressed by the following equations, respectively [47], In order to show normalized horizontal and vertical stresses, Equations (2) and (3) were simplified as follows: On the loading diameter, which is x = 0, Equation (5) is presented in this form: Similarly, on the horizontal diameter, which is y = 0, Equation (6) is presented in this form: In order to show normalized horizontal and vertical stresses, Equations (2) and (3) were simplified as follows: On the loading diameter, which is x = 0, Equation (5) is presented in this form: Similarly, on the horizontal diameter, which is y = 0, Equation (6) is presented in this form: Energies 2020, 13, 1441 7 of 18

Model Assumptions
In Brazilian testing, a cylindrical specimen disc is compressed diametrically (Figure 4a,b). The assumptions [6] are given below: i.
The Brazilian disc is loaded by diametrically opposed uniform pressure at both ends of its diameter over a small strip of the perimeter; ii.
The stress magnitude caused by the internal frictional force between the specimen and loading platens is neglected; iii.
The material is assumed as isotropic, homogeneous, and remains in a linearly elastic zone before causing the brittle failure.

Model Set-up and Inputs
This numerical model ( Figure 5) was built by COMSOL Multiphysics, a multi-field finite element modeling (FEM) software. As it was shown in the Brazilian disc test (BDT), the tensile strength was determined by diametrically compressing the cylindrical rock samples between the loading platens. The basic parameters of the numerical model are given in Table 2.
over a small strip of the perimeter; ii.
The stress magnitude caused by the internal frictional force between the specimen and loading platens is neglected; iii.
The material is assumed as isotropic, homogeneous, and remains in a linearly elastic zone before causing the brittle failure.

Model Set-up and Inputs
This numerical model ( Figure 5) was built by COMSOL Multiphysics, a multi-field finite element modeling (FEM) software. As it was shown in the Brazilian disc test (BDT), the tensile strength was determined by diametrically compressing the cylindrical rock samples between the loading platens. The basic parameters of the numerical model are given in Table 2.
The maximum and minimum mesh sizes were set as 0.8 mm and 0.526 mm, respectively. In COMSOL Multiphysics, we selected the Solid Mechanics Module and Stationary Solver to create the simulation. To obtain the quasi-static stress and strain distribution, three different percentages of peak load (20%, 40%, and 90%) obtained experimentally from Brazilian testing was applied to the top platen.  Figure 5 shows the element and nodal distribution for the Brazilian disc and loading platens. The free quadrilateral mesh was used to define the shapes and boundaries of finite elements and continuous deformation. For this BDT simulation, a load was applied to the rock sample through the upper platen with a constant or variable loading rate over time. Stresses (σx and σy) at any moment can be expressed by Equations (2) and (3), under a known loading P or loading displacement. There are two types of loading platens, which are used very commonly, flat platens and curved platens. In this paper, we applied curved platens in the simulation model, based on our real experimental conditions.  The maximum and minimum mesh sizes were set as 0.8 mm and 0.526 mm, respectively. In COMSOL Multiphysics, we selected the Solid Mechanics Module and Stationary Solver to create the simulation. To obtain the quasi-static stress and strain distribution, three different percentages of peak load (20%, 40%, and 90%) obtained experimentally from Brazilian testing was applied to the top platen. Figure 5 shows the element and nodal distribution for the Brazilian disc and loading platens. The free quadrilateral mesh was used to define the shapes and boundaries of finite elements and continuous deformation. For this BDT simulation, a load was applied to the rock sample through the upper platen with a constant or variable loading rate over time. Stresses (σ x and σ y ) at any moment can be expressed by Equations (2) and (3), under a known loading P or loading displacement. There are two types of loading platens, which are used very commonly, flat platens and curved platens. In this paper, we applied curved platens in the simulation model, based on our real experimental conditions.  Figure 6a shows the dimensionless numerical and analytical stress σ x and σ y on vertical diameter AB during loading. In the center of the rock sample, the normalized stress σ x and σ y were 1 and −3, respectively. Here, tensile stress was defined as positive and compressive stress was negative. In Figure 6a, the maximum value of the numerically horizontal tensile stress, σ x was found at the center of the disc and gradually decreased from the center to the two loading platens, and even became horizontal compressive stress when it came too close to the loading platens.

Numerical vs. Analytical Solutions for Stress Distribution
analytical solutions (Equations (2) and (3)) of BDT, the stress σx and σy will approach negative infinity at the loading area, which did not match the real BDT. To avoid the unrealistic negative infinity in Equations (5) and (6), the value of the minimum dimensionless compressive stress in Figure 6a was not from the loading point which had a distance of a radius from the center of the rock sample, but 95% of the radius instead. In general, the numerical solutions matched the analytical solutions well.
where At is the analytical value, Ft is the numerical simulated value, and n is the number of fitted points. Figure 6b presents the dimensionless numerical and analytical stresses σx and σy on the horizontal diameter CD during loading. The magnitude of the absolute value of the compressive stress, σy was three times the tensile stress σx at the center of the Brazilian disc in the simulation, which matched the analytical solution given by Equations (7) and (8). The magnitude of the absolute value of the stresses σx and σy were 0 at the ends of the horizontal diameter, gradually increasing and getting their maximum values in the center, respectively. Equation (9)  In other words, theoretically, the potential tensile failure area should occur in a region where normalized σx is close to 1. Therefore, the failure stress in the center can be used to compute the tensile strength of the sample. The distribution patterns of σx and σy were symmetric with respect to the yaxis, respectively. The dimensionless stress distribution on the horizontal diameter can be calculated by Equations (5) and (6). The numerical solutions matched the analytical solutions very well in Figure  6b. The MAPEs of σx and σy on the horizontal diameter were 2.95% and 1.58%, respectively. Figure 6c shows the dimensionless stress distribution of σx, σy, and τxy from analytical solutions in the Brazilian disc and Figure 6d shows the respective stresses (stresses both horizontally and vertically, and shear stress) from the numerical solutions. The model input used in the model was based on laboratory testing on Spider Berea sandstone.
From Figures 6c,d, it was evident that the analytical and numerical solution closely matched, while there was a small difference observed near the boundary of the sample. The difference likely On the other hand, the maximum value of the horizontal compressive stress, σ y found at the center of the Brazilian disc gradually decreased from the center to the two loading platens. From the analytical solutions (Equations (2) and (3)) of BDT, the stress σ x and σ y will approach negative infinity at the loading area, which did not match the real BDT. To avoid the unrealistic negative infinity in Equations (5) and (6), the value of the minimum dimensionless compressive stress in Figure 6a was not Energies 2020, 13, 1441 9 of 18 from the loading point which had a distance of a radius from the center of the rock sample, but 95% of the radius instead. In general, the numerical solutions matched the analytical solutions well.
where A t is the analytical value, F t is the numerical simulated value, and n is the number of fitted points. Figure 6b presents the dimensionless numerical and analytical stresses σ x and σ y on the horizontal diameter CD during loading. The magnitude of the absolute value of the compressive stress, σ y was three times the tensile stress σ x at the center of the Brazilian disc in the simulation, which matched the analytical solution given by Equations (7) and (8). The magnitude of the absolute value of the stresses σ x and σ y were 0 at the ends of the horizontal diameter, gradually increasing and getting their maximum values in the center, respectively. Equation (9) is the mean absolute percentage errors (MAPE) formula. The MAPE of σ x and σ y on the vertical diameter were 10.95% and 14.59%, respectively. At the center point, the MAPE of σ x and σ y were 4.02% and 1.12%, respectively.
In other words, theoretically, the potential tensile failure area should occur in a region where normalized σ x is close to 1. Therefore, the failure stress in the center can be used to compute the tensile strength of the sample. The distribution patterns of σ x and σ y were symmetric with respect to the y-axis, respectively. The dimensionless stress distribution on the horizontal diameter can be calculated by Equations (5) and (6). The numerical solutions matched the analytical solutions very well in Figure 6b. The MAPEs of σ x and σ y on the horizontal diameter were 2.95% and 1.58%, respectively. Figure 6c shows the dimensionless stress distribution of σ x , σ y , and τ xy from analytical solutions in the Brazilian disc and Figure 6d shows the respective stresses (stresses both horizontally and vertically, and shear stress) from the numerical solutions. The model input used in the model was based on laboratory testing on Spider Berea sandstone.
From Figure 6c,d, it was evident that the analytical and numerical solution closely matched, while there was a small difference observed near the boundary of the sample. The difference likely occurred because of the stress field; since stresses are theoretical solutions, stress singularity of σ x , σ y , and τ xy (Equations (2)-(4)) was present at loading points that differ from the real experiment. On the other hand, the analytical solutions were always calculated based on exactly applied point load and stress concentration occurred near the loading points whereas the solutions from numerical were always approximate.

Strain Distribution in Homogeneous Rocks from FEM vs. DIC Results
To establish a benchmark, we initiated the experiments with homogeneous Spider Berea sandstone [41] without any lamination, natural fractures, or visible discontinuities. The deposition process of sedimentary rocks is naturally inhomogeneous in nature, but the Spider Berea sandstone analyzed experimentally visibly had no laminations and as such it was assumed as a homogeneous rock.
Images were captured throughout the test until the final failure of the specimen (Figure 7). The strain map created by the image correlation technique aids us to understand the deformation behaviors. The experiment investigated the highest tensile strain in the Spider Berea sandstone sample was along the vertical centerline; with load, the tension on the centerline increased until it was causing a vertical, and central fracture in the peak load of 1277 lb f . behaviors. The experiment investigated the highest tensile strain in the Spider Berea sandstone sample was along the vertical centerline; with load, the tension on the centerline increased until it was causing a vertical, and central fracture in the peak load of 1277 lbf. A few snapshots at different percentage peak loads (20%, 40%, and 90%) during testing are illustrated here. The peak load for Berea sandstone was obtained from the experimental investigation conducted using the Instron load frame (Figure 3). The horizontal (εxx), vertical (εyy) and shear strain (εxy) is presented at these three defined peak loads (Figure 7), which is capable of explaining the mechanism involved in deformation and failure including tension, compression or shear.
The peak obtained from the experiment was 1277 lbf and three different peak loads (255, 510 and 1150 lbf respectively) were selected based on the elastic limit. Horizontal strain (εxx) distribution for three different peak loads obtained from DIC were compared with simulation results generated by COMSOL Multiphysics, which is presented in Figure 8.
From Figure 8, good agreement was demonstrated between horizontal strain distributions from both DIC and the model. The simulation results showed that compressive strain was generated near the loading point or boundary, which was relatively less visible in DIC as it could not compute deformation up to half subset size near the boundary of the speckled region. The results also substantiated there was a reasonably good match obtained within the elastic response (20% and 40% of peak load). As the load increased and the sample reached closer to failure, the deviation in deformation along with the centerline increased as shown for 90% of the peak load.
The DIC color maps displayed some sprinkle nature in Figure 8. At 20% of the peak load, it was evident that there were some dark-yellow (black circle) and light-yellow (black rectangle) features that looked like noise in the measurements during loading. According to the standard, DIC noises were within 0.01% strain. Dark-yellow dots were above these noise thresholds and represent a real physical deformation that was measured using DIC. The light-yellow features have low strain values and are very difficult to measure using DIC as well as conventional techniques. The light-yellow features may be the noise or deformations that cannot be captured using DIC with high fidelity. A few snapshots at different percentage peak loads (20%, 40%, and 90%) during testing are illustrated here. The peak load for Berea sandstone was obtained from the experimental investigation conducted using the Instron load frame (Figure 3). The horizontal (ε xx ), vertical (ε yy ) and shear strain (ε xy ) is presented at these three defined peak loads (Figure 7), which is capable of explaining the mechanism involved in deformation and failure including tension, compression or shear.
The peak obtained from the experiment was 1277 lb f and three different peak loads (255, 510 and 1150 lb f respectively) were selected based on the elastic limit. Horizontal strain (ε xx ) distribution for three different peak loads obtained from DIC were compared with simulation results generated by COMSOL Multiphysics, which is presented in Figure 8. At 40% of the peak load, most yellow dots in the DIC map had a sprinkled look and were not as smooth as in the simulation. This could probably relate to the porosity variation and nonhomogeneity of the Spider Berea rock.
Moreover, at 90% of the peak load, the DIC map was similar to the simulation magnitude. It substantiated that the homogeneous simulation could capture strain deformation events of rock under diametrical compression.
Similarly, vertical strain (εyy) distribution for three different peak loads also obtained from DIC was compared with the simulation results, which were presented in Figure 9. The strain range was selected to get the maximum visibility of the deformation behaviors. The vertical strain (εyy) distribution also matched close to the DIC generated value. From Figure 8, good agreement was demonstrated between horizontal strain distributions from both DIC and the model. The simulation results showed that compressive strain was generated near the loading point or boundary, which was relatively less visible in DIC as it could not compute deformation up to half subset size near the boundary of the speckled region. The results also substantiated there was a reasonably good match obtained within the elastic response (20% and 40% of peak load). As the load increased and the sample reached closer to failure, the deviation in deformation along with the centerline increased as shown for 90% of the peak load.
The DIC color maps displayed some sprinkle nature in Figure 8. At 20% of the peak load, it was evident that there were some dark-yellow (black circle) and light-yellow (black rectangle) features that looked like noise in the measurements during loading. According to the standard, DIC noises were within 0.01% strain. Dark-yellow dots were above these noise thresholds and represent a real physical deformation that was measured using DIC. The light-yellow features have low strain values and are very difficult to measure using DIC as well as conventional techniques. The light-yellow features may be the noise or deformations that cannot be captured using DIC with high fidelity.
At 40% of the peak load, most yellow dots in the DIC map had a sprinkled look and were not as smooth as in the simulation. This could probably relate to the porosity variation and non-homogeneity of the Spider Berea rock.
Moreover, at 90% of the peak load, the DIC map was similar to the simulation magnitude. It substantiated that the homogeneous simulation could capture strain deformation events of rock under diametrical compression.
Similarly, vertical strain (ε yy ) distribution for three different peak loads also obtained from DIC was compared with the simulation results, which were presented in Figure 9. The strain range was selected to get the maximum visibility of the deformation behaviors. The vertical strain (ε yy ) distribution also matched close to the DIC generated value. At 40% of the peak load, most yellow dots in the DIC map had a sprinkled look and were not as smooth as in the simulation. This could probably relate to the porosity variation and nonhomogeneity of the Spider Berea rock.
Moreover, at 90% of the peak load, the DIC map was similar to the simulation magnitude. It substantiated that the homogeneous simulation could capture strain deformation events of rock under diametrical compression.
Similarly, vertical strain (εyy) distribution for three different peak loads also obtained from DIC was compared with the simulation results, which were presented in Figure 9. The strain range was selected to get the maximum visibility of the deformation behaviors. The vertical strain (εyy) distribution also matched close to the DIC generated value. Both the DIC and FEM deformation patterns were very similar for 20% of the peak load in the vertical strain map shown in Figure 9. At 40% and 90% of the peak load, the strain maps showed some blue horizontal irregular lines, which probably occurred due to the compression of the soft Both the DIC and FEM deformation patterns were very similar for 20% of the peak load in the vertical strain map shown in Figure 9. At 40% and 90% of the peak load, the strain maps showed some blue horizontal irregular lines, which probably occurred due to the compression of the soft lamination inherently present in the rock. The deformation of 90% of the peak load reasonably matched with the homogeneous simulation.
Likewise, the shear strain (ε xy ) distribution for the aforementioned three different peak loads obtained from DIC was also compared with the COMSOL model and is illustrated in Figure 10. The results have shown that the shear strain was less than the axial strain, which substantiates that the shear strain was not responsible for the failure of isotropic and homogeneous rock. The shear strain (ε xy ) distribution also matched within 2% of the DIC generated value. Figure 10 also depicted the sprinkle nature on the color map, which most probably occurred owing to the porosity, soft lamination, and non-homogeneity nature of the sample, but the overall deformation matched closely. There were some other reasons such as the vibration of the load frame, temperature from the heating of the camera, camera distortion, etc. which may have caused noise that resulted in the deviation in the results during the experiment [48] but these can be minimized by Energies 2020, 13, 1441 12 of 18 calibration. The above DIC and simulation results revealed that the homogeneous assumption in the simulation can be acceptable.
obtained from DIC was also compared with the COMSOL model and is illustrated in Figure 10. The results have shown that the shear strain was less than the axial strain, which substantiates that the shear strain was not responsible for the failure of isotropic and homogeneous rock. The shear strain (εxy) distribution also matched within 2% of the DIC generated value. Figure 10 also depicted the sprinkle nature on the color map, which most probably occurred owing to the porosity, soft lamination, and non-homogeneity nature of the sample, but the overall deformation matched closely. There were some other reasons such as the vibration of the load frame, temperature from the heating of the camera, camera distortion, etc. which may have caused noise that resulted in the deviation in the results during the experiment [48] but these can be minimized by calibration. The above DIC and simulation results revealed that the homogeneous assumption in the simulation can be acceptable.

Vertical Displacement from Numerical Model vs. DIC Results
At three different portions of applied peak loads, the changes in displacement along the central vertical axis during loading were also investigated to understand the deformation behaviors and compared with results from the numerical model which were presented in terms of displacement contour in Figure 11. Figure 12 demonstrated that 20% and 40% of the peak values showed very good agreement both in numerical and DIC generated displacement magnitude. For 90% of the peak load, the MAPE was investigated within 1.21% from the experimental value for vertical displacement using Equation (9). Overall, vertical displacement was investigated as repetitive and close to the experimental value.

Vertical Displacement from Numerical Model vs. DIC Results
At three different portions of applied peak loads, the changes in displacement along the central vertical axis during loading were also investigated to understand the deformation behaviors and compared with results from the numerical model which were presented in terms of displacement contour in Figure 11.   Figure 12 demonstrated that 20% and 40% of the peak values showed very good agreement both in numerical and DIC generated displacement magnitude. For 90% of the peak load, the MAPE was investigated within 1.21% from the experimental value for vertical displacement using Equation (9). Overall, vertical displacement was investigated as repetitive and close to the experimental value.

348
A similar investigation was conducted on the changes in horizontal displacement and a very 349 reasonable agreement in both approaches was found ( Figure 13).

Horizontal Displacements from the Numerical Model vs. DIC Results
A similar investigation was conducted on the changes in horizontal displacement and a very reasonable agreement in both approaches was found ( Figure 13).
Energies 2020, 13, x FOR PEER REVIEW 14 of 20 A similar investigation was conducted on the changes in horizontal displacement and a very reasonable agreement in both approaches was found ( Figure 13).   Table 3); 90% of the final peak values showed a comparatively higher difference. Overall, horizontal displacement was investigated as repetitive and closer to the experimental value.   Table 3); 90% of the final peak values showed a comparatively higher difference. Overall, horizontal displacement was investigated as repetitive and closer to the experimental value. Figure 14 demonstrated that 20% and 40% of the final peak values showed an excellent agreement both in numerical simulation and DIC results. The deviation increased as the load increased ( Table 3); 90% of the final peak values showed a comparatively higher difference. Overall, horizontal displacement was investigated as repetitive and closer to the experimental value.  From the R 2 value, it was evident that the deviation along the vertical centerline (Dy) was bigger than that of the horizontal (Dx). This is probably because of the existence of non-linearity in the specimen, in particular, because of soft lamination and the bedding nature. This linear elastic numerical model did not take into account such a non-linear effect.
Considering the comparison between two distinct stages of the failure process, such as, right before failure and during failure of the sample, can be utilized to investigate the possible crack initiation locations as shown in Figure 15. In the images right before the failure stage, it was evident that the intensity of horizontal strain accumulation (ε xx ) was greater in the central region compared to the boundary of the sample. During failure, central crack propagation was observed up to the region where strain accumulated greater, and not near the loading points. Similarly, the horizontal displacement (Dx) plot right before failure showed that the central region was being pulled apart (by tension) in two opposite directions, while the horizontal displacement was less in the contact regions. Therefore, the generated cracks that were initiated from the center or central region of the sample substantiated the validity of the indirect tensile testing based on the fracture mechanics theory. Comparison of analytical, numerical and experimental results using DIC was in considerable agreement, so the simulation model based on the experimental investigation is effective at predicting the deformation behavior of rocks. The wider difference between the two approaches was found near to the loading points, as it was this zone that was more prone to stress concentrations, which could be addressed by future research. Energies 2020, 13, x FOR PEER REVIEW 16 of 20 The limitation of Brazilian experiments is that it cannot imitate the real subsurface condition which is fluid-driven, but it aids to get an improved understanding of induced fracture geometry when the formation is geologically complex with different discontinuities. Where the complex rock failure process is critical to understand by experimental investigation, the numerical model validated by experimental properties better characterizes the deformation properties of the rock. Overall, this study presented a benchmark for characterizing the deformation attributes for homogeneous rocks and the variation due to the non-homogeneity of rock, will be captured by our future research.

Conclusions
In this study, numerical simulations of Brazilian testing based on the finite element method using COMSOL Multiphysics was conducted on the Berea sandstone and fracture behaviors were compared with experimental investigation using the DIC. For analyzing the deformation behaviors, the following conclusions were made: i.
The highest tensile strain investigated along the vertical centerline in the Spider Berea sandstone sample; with load, the tension on the centerline increases until it is causing a vertical, and central fracture; ii.
Linear elastic quasi-static numerical results of Brazilian testing reasonably match the analytical solutions. The MAPEs of σx and σy on the vertical diameter are 10.95% and 14.59%, respectively; iii.
The horizontal (εxx), vertical (εyy) and shear strain (εxy) were investigated for 20%, 40% and 90% of the peak load obtained from the experimental investigation. The simulation results obtained from the numerical model matched reasonably with DIC generated strain maps except for the proximity of the boundary as DIC maps are in limitations encountering the boundary effect. Therefore, a validated simulation model can be used to predict the boundary effect in Brazilian testing with a different condition for which experimental validation is critical; iv.
Changes in vertical (Dy) and horizontal (Dx) displacement in the along centerline was investigated. The numerical results obtained were repetitive and closely match the experimental The limitation of Brazilian experiments is that it cannot imitate the real subsurface condition which is fluid-driven, but it aids to get an improved understanding of induced fracture geometry when the formation is geologically complex with different discontinuities. Where the complex rock failure process is critical to understand by experimental investigation, the numerical model validated by experimental properties better characterizes the deformation properties of the rock. Overall, this study presented a benchmark for characterizing the deformation attributes for homogeneous rocks and the variation due to the non-homogeneity of rock, will be captured by our future research.

Conclusions
In this study, numerical simulations of Brazilian testing based on the finite element method using COMSOL Multiphysics was conducted on the Berea sandstone and fracture behaviors were compared with experimental investigation using the DIC. For analyzing the deformation behaviors, the following conclusions were made: i.
The highest tensile strain investigated along the vertical centerline in the Spider Berea sandstone sample; with load, the tension on the centerline increases until it is causing a vertical, and central fracture; ii.
Linear elastic quasi-static numerical results of Brazilian testing reasonably match the analytical solutions. The MAPEs of σ x and σ y on the vertical diameter are 10.95% and 14.59%, respectively; iii.
The horizontal (ε xx ), vertical (ε yy ) and shear strain (ε xy ) were investigated for 20%, 40% and 90% of the peak load obtained from the experimental investigation. The simulation results obtained from the numerical model matched reasonably with DIC generated strain maps except for the proximity of the boundary as DIC maps are in limitations encountering the boundary effect. Therefore, a validated simulation model can be used to predict the boundary effect in Brazilian testing with a different condition for which experimental validation is critical; iv. Changes in vertical (Dy) and horizontal (Dx) displacement in the along centerline was investigated. The numerical results obtained were repetitive and closely match the experimental DIC value. Dx comparison was better matching than Dy due to non-linearity arising from the bedding plane compression, which could not be captured by the linear elastic model; v.
Simulation results of 20% and 40% of the peak loads for vertical (Dy) and horizontal (Dx) displacement are in very good agreement till the load increased near to the failure (90% of the peak load) the DIC value deviates more from simulation value. These observations revealed that changes in both vertical and horizontal displacement are very little in the elastic limit. vi.
DIC technique is capable of producing a strain map across the rock surface and can capture the effects of porosity variation, soft lamination, and other non-homogeneity.

Conflicts of Interest:
There is no conflict of interest in our research paper. Stress in the x-direction σ y Stress in the y-direction τ xy Shear stress ε xx Strain in the x-direction ε yy Strain in the y-direction ε xy Shear strain Dx

Abbreviations
Horizontal strain along the centerline Dy Vertical strain along the centerline