Four-point bending test on a middle strength rock : numerical and experimental investigations

Developing a reliable numerical modelling technique is considered as challenge for fracture assessment of the geological materials, which are much subjected to hydrostatic pressure. For this purpose, the mechanical behaviour and the fracture pattern of a middle strength rock material, called Pietra Serena sandstone, is investigated both numerically and experimentally under a Four-Point Bending (also called Flexural) testing program. For the numerical approach, an innovative method, namely FEMcoupled to-SPH, is exploited due to its capabilities in dealing with rock mechanics related issues. Two different material models, which are the Karagozian and Case Concrete (KCC) and the Extended (Linear) DruckerPrager, are exerted to assess their capabilities. The Flexural strength and the crack initiation area are studied based on the state of the stress in various parts of the specimen in both models, and finally the results obtained from the numerical models are compared with the data obtained from the experimental tests in order to assess the capability of the modelling approach.


INTRODUCTION
ock materials are highly dependent on the hydrostatic loading pressure, moreover their tensile strength is considerably weaker than the compressive one.Thus, the characterization of the tensile behaviour is of great importance in geophysical applications.The tensile behaviour of rock materials, which is also one of the most important parameters to predict the rock's boreability, can be determined in several ways, i.e. by a direct tensile test, a Brazilian splitting disk test, a (three-or four-point) bending test, etc [1].The direct tensile test seems to be the standard way to determine the maximum principal tensile strength of rock materials.However, due to several issues, including the complicated set up, the non-uniformity and the perturbation in the uniaxial stress field introduced by specimen grips and/or slight imperfections in the sample preparation and the material inhomogeneity, etc., its usage is inconvenient [2][3][4].Therefore, rock engineers perform the Brazilian and the Flexural (four-point bending) test to investigate the tensile behaviour.These tests are designed to investigate indirectly a normal tensile stress state (at least) in a specific portion of the specimen.The numerical and experimental investigation of the mechanical behaviour of a middle strength rock material which is subjected to a Four-Point Bending (or Flexural) test is the aim of this research study.The Four-Point Bending test is an indirect way to estimate the tensile strength, which consists of a beam in flexure.The strength of the beam, in this test, is expressed mainly in terms of modulus of rupture, which usually tends to overestimate the tensile strength up to one-hundred percent [4].This overestimation is mainly caused by the hypothesis of the test which R assumes that the mechanical characteristics, i.e. the stress-strain behaviour, are linear all the way through the critical crosssectional area of the beam [5].Due to the sensitivity of this test to the boundary conditions, it may be inaccurate in large displacement (which is similar to the direct tensile test) [6].In addition to these difficulties, the strength of quasi-brittle materials, e.g.rock, is considered to be dependent on their size and scale [7][8][9].However, the Flexural test is considered as the simplest method to investigate the creep (time-dependent behaviour) of rock [10], and accordingly, represents an interesting alternative for the investigation of the stress-strain relationships [11,12].The experimental campaign designed for this research study follows the protocols of the ASTM standard and the results can therefore be easily compared with similar studies, i.e. on other materials.The mechanical properties of several types of rock materials, including sandstones have been studied in a large variety of research.Among them, the Berea sandstone which is deposited sub-aqueously as offshore beds in a considerable number of well drilling applications, is widely investigated by the oil and gas industries [13,14].It is a sedimentary rock with high porosity and mostly composed of sand-sized grains [15], however, it is not a highly accessible material resulting in a relatively high expense for testing purposes.An extensive literature review [16][17][18] exposed the presence of another rock material, which is called Pietra Serena sandstone, with similar mechanical characteristics to Berea sandstone.These similar properties led the authors to perform the experimental tests on the Pietra Serena sandstone.It is anticipated that the experimental data of the flexural test on the Pietra Serena can be used to have an outlook of the behaviour of Berea sandstone.The development of a reliable numerical method in conjunction with an accurate constitutive material model has been considered an essential tool in the stress analyses.Several numerical modelling techniques have been recently implemented [19][20][21][22][23] to simulate the rock materials.One of the most common and accurate numerical simulation techniques, which is implemented in research fields as well as in industrial applications, is the Lagrangian Finite Element Method (FEM).However, this method cannot appropriately deal with large deformations and tearing, which are often present in the numerical modelling of fractured rock.The Smooth Particles Hydrodynamics (SPH), on the other hand, is a mesh-free Lagrangian method that discretizes a system into a number of "mesh-points" (or particles) carrying the field variables.Due to the not fixed nodal connectivity of this method, the SPH is able to cope with highly distorted elements [21,24,25].The FEM is however more efficient in terms of accuracy and computation time when SPH particles are used.Therefore, an advanced technique is exploited, inspired by the research study of the same authors of the present paper, Bresciani et al., [25] to take advantage of both the Lagrangian FEM (before the occurrence of high distortion) and the SPH methods to deal with large deformation, mesh distortion, etc.In this numerical model, which is called FEM-coupled to-SPH, the specimen is initially modelled by Lagrangian 3D elements and subsequently by means of an eroding algorithm.The elements which reach a specific failure level are eroded, and subsequently these eroded elements are transformed to SPH particles with the same mechanical properties.Two commercial numerical solvers, LS-DYNA and ABAQUS, are used to replicate the experimental tests in conjunction with two constitutive material models: the Karagozian and Case Concrete (KCC) model and the Extended (Linear) Drucker-Prager (LDP) model.The KCC is an advanced material model, developed by Malvar et al. [26][27][28][29], available in LS-DYNA, that decouples the volumetric and deviatoric responses.This material model consists of three-independent failure surfaces to determine the accumulated damage.The LDP, on the other hand, is based on the conventional Drucker-Prager model.This material model which is available in ABAQUS, takes advantage of some improvements, i.e. the flow rule plasticity [30] and is particularly interesting due to the potential use of different shapes of the yield function (linear, hyperbolic and exponential), and additional cap yield function and the possibility to be used in conjunction with an equation of state.However, due to the lack of a tension (or pressure) cut-off level, the numerical results have to be critically evaluated.The numerical results of the KCC and LDP are thus compared in detail to show the reliability of the models.The article is divided into the following sections.The methods and the results obtained during the experimental tests are reported in section 2. In section 3, the theories of the KCC and the LDP material models are discussed in detail.Then the experimental configuration suggested by the ASTM is replicated in LS-DYNA and ABAQUS by the FEM-coupled to-SPH method in conjunction with the KCC and the LDP material models, respectively.The numerical results of these models are then compared with each other and discussed in section 4 by the further comparison of the numerical modelling results with the experimental data.

EXPERIMENTAL TEST
he experimental configuration for the Flexural test is designed based on the protocols of the ASTM standard [31].This configuration consists of a rectangular cubic specimen which is supported by two fixed rollers near the end of its length span (see Fig. 1a).Thus, the specimen is loaded vertically by means of two compressive rollers at a certain T distance from the center of the specimen.This symmetrical configuration causes nominally zero shear forces, and accordingly constant bending moment between the two compressive rollers.The normal tensile and compressive stresses appear at the top and bottom of this middle span, respectively.According to the beam theory, the maximum principal stress corresponding to the ultimate loading value can be determined, which is called the flexural strength and gives a rough approximation of the principal tensile strength.The flexural strength tends to overestimate the tensile strength, because the measuring process considers a linear relationship (as the stress-strain behaviour of the material) at the beam critical cross-section.Furthermore, all the materials have a certain amount of anisotropic level in their structure [4].The flexural strength .flex


[MPa], which is given by the Eq. ( 1), can be considered as a parameter to validate the numerical models.
where, W was measured as the maximum applied force.The experimental tests within this study were performed on five specimens of Pietra Serena sandstone with the same geometries.The span length L, width b and height d of all the specimens were equal to 318, 102 and 32 [mm], respectively.However, the total length of the specimen was measured as 381 [mm].Two pairs of steel rods were embedded to the testing apparatus (see Fig. 1).The axes of lower rods were fixed to the bed of testing machine while the upper rods were displacement controlled by means of a compressive platen.According to [31] the speed was set to 0.2 [mm/min] in order to apply the load at a uniform stress rate of 4.14 [MPa/min].This study is aimed to investigate the mechanical response of the rock based on both the maximum loading and the displacement.The applied load of the compressive platen was measured automatically by the load cell of the apparatus.It is not mentioned how to measure the displacement of the specimen on the ASTM standard.It was decided to measure the displacement of the specimen by means of a displacement gage contacting extensometer.As can be seen in Fig. 2, the flexible tip of the extensometer is located at the center of the lower face of specimen.However, since the extensometer, itself, is fixed to the fixture of the testing machine, the displacements of all the components in between are measured.Therefore, in the configuration expressed in Fig. 2, the experimental data provided by the extensometer contains the displacements of the specimen, the fixed rollers and the steel blocks between the rollers and the fixture of the apparatus.This set of data therefore doesn't represent the actual displacement of the specimen, but since the mechanical properties of the rollers and the steel blocks are known, they can be replicated in the numerical simulations as well.This set up is thus a convenient way to record the experimental data of the flexural test in terms of displacement.The broken specimens after performing the flexural test are indicated in Fig. 3.In all of the specimens, cracks initiate at the lower part, which is subjected to tension stresses and then propagate in an upward manner through the depth.Also, the cracks are located under (and close to) the section of the specimens which are in contact with the moving rods.The maximum load, displacement and the flexural strength are reported in Tab. 1, separately for each specimen.Average values, 95% Confidence Interval ("95% CI") of the average value and standard deviations has been also calculated from experimental data.CI is an interval that try to estimate range of the average value of the population from the sample data.Due to the large variability of the mechanical behavior of rock this range is more representative than the single average value obtained from the 5 experimental tests.Fig. 4

NUMERICAL MODELLING
any numerical modelling techniques have been developed for stress analyses of solid mechanics in the continuum domain.Among these, the non-linear Lagrangian Finite Element Method has its own privileges for a number of reasons, i.e. the accuracy and the convenient time consumption cost.However, a disadvantage lies within the reduced performance of grid-based numerical methods to deal with highly distorted elements and fragmentation.At present, the Smooth Particle Hydrodynamics (SPH) is one of the most convenient numerical methods to cope with problems related to large deformations and fracture.This method is a numerical mesh-free approach that represents the state of a system by a set of particles that hold the properties of the material.The particles interact in a range which is controlled by a kernel (smoothing) function, so that this kernel approximation for any two given particles, i and j, can be expressed by Eq. ( 2).
where, the W in Eq. ( 2) is the smooth function which can be determined by Eq. (3).
The i x and j x are the coordinates of the particles in the problem domain Ω and h is the distance between the two particles.The parameters d and θ in Eq. ( 2) represent the number of the space dimension and an auxiliary function, respectively (see

Load [kN]
Specimen K1 Specimen K2 Specimen K3 Specimen K4 Specimen K5 M Fig. 5).Because the connectivity between these particles is considered as a part of the computation procedure, a straightforward handling of problems with large deformations is allowed by using the SPH method [33].A numerical technique, namely FEM-coupled to-SPH, was developed, inspired by [25], to take advantage of both the FEM and the SPH methods.The Lagrangian meshed elements were eroded after reaching a certain criterion by using this technique, and they were adaptively transformed to SPH particles that inherit identical mechanical properties with the failed elements.The implementation of this technique in LS-DYNA was carried out by calling the keyword called ADAPTIVE_SOLID_TO_SPH.However, since many of the material keywords available in the LS-DYNA material library, i.e. the ones used in this study (KCC), do not have any internal eroding algorithm, an external eroding algorithm is required to be implemented in conjunction with this keyword.The element erosion of the models of this article was obtained by using MAT_ADD_EROSION and specifying a certain value for the maximum shear strain (EPSSH).Therefore, each hexahedral solid element who meets this criterion will be eroded, and then by defining the two input parameters of ADAPTIVE_SOLID_TO_SPH keyword as ICPL=1 and IOPT=1, the software automatically replaces those eroded elements with a certain number of SPH particles.The number of particles can be controlled by the users, i.e. it can be 1, 8 or 27 in case of hexahedral elements (see Fig. 6).Within this study, the failed solid elements were converted into one SPH particle to keep the time consumption cost low.This method can be implemented in ABAQUS by choosing the option "conversion to particles" from "element type" panel at the mesh module.Unlike the LS-DYNA, an external eroding algorithm in never required since the elements which meet a user-specified criterion are automatically transformed to a certain number (from 1 up to 343) of SPH particles.The maximum principle strain is considered as the conversion criterion for all the models developed by ABAQUS in this study.

Karagozian and Case Concrete (KCC) model
The studies on the nonlinear concrete material models that were implemented by the Lagrangian finite element code DYNA3D at 1992 (reported in [34]), showed that the material model  ( , , ) , where 1 I is the first principal invariant of the Cauchy stress, 2 J and 3 J are the second and third principal invariants of the deviatoric part of the Cauchy stress, respectively.It is more convenient to use a cylindrical set of coordinate system for the cohesive frictional material since the principal stresses can be immediately recast into this format as Eq. ( 4).The formulation of KCC model is explained by means of this cylindrical coordinate system as a clearer way of guidance for readers.
where, the  ,  and  are the Haigh-Westergaard coordinates defined as Eq. ( 5).
  The failure function of the KCC material model can be derived by Eq. ( 6) in terms of the Haigh-Westergaard stress invariants.
where, ( , ) c f   corresponds with the failure surface in the triaxial compression state of stress, when the Lode angle θ is equal to 60  , and the non-dimensional function ˆ[ ( ), ] r p   is the ratio between the current radius of the failure surface and the compressive meridian.The Eq. ( 6) is discussed in detail below.The Karagozian and Case concrete model takes advantage of the three-fixed independent failure surfaces in the compressive meridian (    plane), which correspond to the yield, the ultimate and the residual strength of the material.Hence, the function ( , ) c f   is defined for each one of these failure surfaces separately by Eqs. ( 7), ( 8) and (9).
The stress invariants for an axisymmetric loading condition, when 1 a    is the axial stress and 2 are the lateral stresses, can be re-written as Eq.(10).
Therefore, the equations of the three-failure surfaces for the axisymmetric state of stress can be written as the Eqs.( 11), ( 12) and ( 13).
where, the i a -parameters are the user-defined input parameters to define the failure surfaces in the compressive meridian.These parameters can be calibrated based on the experimental data of the triaxial compression test by means of a curvefitting approach, e.g. a genetic algorithm.The parameter p is the pressure (positive in compression) or the mean stress equal to 1  As written in the Eq. ( 6), the KCC material model considers the effect of the third invariant, i.e. the Lode angle θ, by means of the function ˆ[ ( ), ] r p   .This function was derived from the Eq. ( 14), which is the shape of the failure criterion in the deviatoric plane, proposed by Willam and Warnke [35].
where, ( ) r  determines the distance of the failure surface at the deviatoric section by considering the effect of the Lode angle θ.The parameters c r and t r express the distances of failure surfaces from the hydrostatic axis at the compressive and tensile meridian, respectively.The deviatoric plane of a Willam-Warnke failure model is indicated in Fig. 8.In this figure, the TXC, TXE and SHR stand for triaxial compression, triaxial tension and pure shear respectively.The ˆ[ ( ), ] r p   , which is the ratio between the current radius of the failure surface ( ) r  and the distance of the failure surfaces from the hydrostatic axis at the compressive meridian c r , is computed by means of the Eq. ( 15).This equation was obtained by dividing both sides of Eq. ( 14) by c r .In order to present the term ( ) p  , which is a strength index of brittle material related to the confining pressure that a material is subjected to it and equal to The fact that r is just a function of ( ) p  and θ, and the lode angle can be determined based on the loading conditions, implies the role of ( ) p  for computational purposes.It means that the implementation of the three invariant failure surfaces is completed by means of this parameter ( ) p  .This parameter generally depends on the hydrostatic stress and can be obtained empirically.Malvar et al. in [28] defined this parameter as a linear piecewise function on the full range of pressure according to Eq. ( 16).
Where c f  is the unconfined compression strength, t f is the principal tensile strength and α is an experimental parameter related to the biaxial compression test.According to the Eq. ( 16), ( ) p  varies from 1⁄2 to 1, which is in accordance with the experimental data previously obtained.It also indicates that 8.45 c p f   is the transition point in which the compression meridian is equal to tension one, and accordingly from this point onwards, there is a circular failure surface on the deviatoric plan section.Moreover, it considers a value equal to 1⁄2 for the negative range of pressures.It is worth mentioning that this function was implemented in LS-DYNA and no input is required of the users.One of the most noteworthy features of the KCC model is decoupling the shear and compaction behaviour of the materials, which means that this model treats the deviatoric and volumetric responses separately.The deviatoric response is characterized by the migration of the current stress state between the fixed failure surfaces, while the response to pressure is defined by an equation of state as a function of volumetric strain increments.The damage accumulation of the KCC model, and accordingly the current failure surfaces are expressed schematically in Fig. 9.As can be seen in Fig. 9a, the state of stress is determined by a linear interpolation between the three failure surfaces.The stress-strain diagram corresponding to the unconfined compression test is indicated in Fig. 9c, which is determined by association of a damage accumulation function (shown in Fig. 9b.The response of the material to the initial loading, phase I in Fig. 9c, is considered as a linear elastic deformation before reaching point 1 in Fig. 9a.The current failure surface . .

( )
Based on the level of confining pressure, a softening response occurs after reaching the maximum strength and before obtaining a residual strength.The current state of stress corresponds to this area (between point 2 and point 3) is derived according to Eq. ( 18).

( )
where η is the KCC damage parameter used to determine the relative location of the current failure surface, which is a function of another parameter, called the accumulated effective plastic strain λ.As can be seen in Fig. 9b, the damage function is described so that initially and prior to the occurrence of any plasticity responses, the value of η is set to zero.This material model implements shear damage accumulation by including the rate effect enhancement and different scenarios in tension and compression as Eq. ( 19). where, is the deviatoric part of strain).The parameter s  expresses the strain rate effect user-defined factor; if it is equal to 0, the strain-rate effects is toggled off, while if it is equal to 100, the strain-rate effects is included.The parameter f r is the dynamic increase factor that accounts for the strain rate effects.Parameters 1 b and 2 b govern the softening in compression and uniaxial tensile strain, respectively, whereas, 3 b , affects the triaxial tensile strain softening.These b-parameters can be determined by iteration until the value of the fracture energy, f G , converges for a specified characteristic length, which is associated with the localization width (i.e. the width of the localization path transverse to the crack advance).
To describe the rock compaction behaviour in LS-DYNA, this *MAT_072R3 is used in conjunction with an Equation-of-State (EOS), namely *EOS_TABULATED_COMPACTION.This EOS provides a function of the current and previous volumetric strain v  for the current element pressure p.The stress tensor can be computed as being a point on a movable surface by means of known pressure.This surface can be a yield or a failure surface.This function should be specified by users as a series of , , v p K  sets, where K is the bulk modulus correspondent to the different  pairs [20].These series can be generated internally by the automatic input generation mode of this material model.In case of an essential change of stiffness of a material, the linear bulk modulus which is defined as the 20, can be adjusted.
The most significant improvement provided by the third release of this material model is the automatic input parameter generation capability, based solely on the unconfined compression strength [29].This makes the use of the KCC model easier and accessible to most users, since there is no need to carry out extensive material characterization tests.The KCC model was originally developed to analyze the mechanical behaviour of concrete, and although the response of sandstone is expected to be similar to concrete, the results of the automatic input generation mode are only a rough estimation for sandstones.However, this automatic method was used initially within this study due to the overwhelming number of input parameters which are difficult to obtain from current resources.Apart from the numerical results which are obtained directly from the automatic mode, LS-DYNA also generates all the input parameters automatically and writes them into the "MESSAG" file.Therefore, some of these parameters (i.e. the b-parameters) can be used in addition to the other parameters provided for the users by other resources.
Previous studies [36,37] showed that the mechanical behaviour of Berea sandstone is expected to be very similar to the Pietra Serena.Since the experimental data related to the triaxial compression test of Pietra Serena sandstone is not available (which is required for calibrating the i a -parameters), it was decided to use the experimental data of the Berea sandstone for the calibration the KCC material model.The initial input parameters used for the automatic mode were obtained from [38] and reported in Tab. 3.

Density, RO [ton/mm 3 ]
Poisson ratio, PR UCS, A0  The parameters RSIZE and UCF are unit conversion factors and the NOUT is called the "output selector for effective plastic strain".According to [39], when NOUT=2, the quantity labelled as "plastic strain" by LS-PrePost is actually the quantity that describes the "scaled damage measure, δ" which varies from zero to two.When the amount of δ is still lower than one, the elements of the part modelled by the KCC, don't reach the yield limit.When this amount is equal to two, the corresponding elements meet the residual failure level.Only at the automatic mode of MAT_072R3 keyword (when a maximum of six parameters are defined) the A0 parameter must be defined as a negative number to represent the unconfined compression strength.For instance, the -62 in Tab.Where, the parameters ω, Edrop and LOC-WIDTH are the frictional dilatancy, post peak dilatancy and three times of the maximum aggregate diameter, respectively.Almost all of the EOS tabular data obtained from the automatic input generation mode, were also used for the full input calibration mode.Only the "Pressure02" parameter from EOS keyword was changed to 26.1 [MPa] to reach the same bulk modulus (and accordingly the elastic behaviour) as the Berea sandstone.The tensile strength parameter was investigated in another recent study by the same authors [37] where the Brazilian tensile test was performed on Pietra Serena and the tensile strength was determined as 5.81 [MPa].
The pairs of η -λ parameters were changed to the values reported in [40].A new set of i a -parameters were computed by the least square curve fitting method from the experimental triaxial compression tests [20].These experimental data which were performed at four different confining pressure levels are reported in [41].The i a -parameters are reported in Tab. 5 (note that the A0 parameter here is a positive value).

Extended (Linear) Drucker-Prager model
This is an advanced material model implemented by ABAQUS, which is based on the conventional Drucker-Prager constitutive law.There are three forms of the extended Drucker-Prager material model developed in ABAQUS which are the linear, the hyperbolic and the general exponential forms [30].Within this research work, the linear form of this model was investigated.The linear Drucker-Prager (LDP) model consists of a pressure-dependent three-invariant failure surface which obtains the opportunity of having the non-circular yield surface (unlike the conventional Drucker-Prager) in the deviatoric plane.This model also takes advantage of the associated inelastic flow and hardening (or softening) rules, and different dilation and friction angles.The general form of the failure surface in the Haigh-Westergaard coordinate system is given by Eq. ( 21).tan ( , , ) ( , ) ( ) The effect of the Lode angle θ is considered by ( , ) t   according to the Eq. ( 22).
where the material parameter K is defined as the ratio of yield stresses at triaxial tension to triaxial compression, which means it controls the shape of the linear Drucker-Prager yield function at the deviatoric plane.As can be seen in Fig. 10b, if K=1.0, a circular failure function presents in the deviatoric plane, which is the same as the conventional Drucker-Prager.
It is proved that in order to verify the convexity of the yield function, the value of parameter K should be in the range of 0.778 to 1.0.By means of an approach to match the parameters of the Mohr-Coulomb and linear Drucker-Prager (for materials with low friction angles), K can be defined according to Eq. ( 23).
where, φ is the internal friction angle used in Mohr-Coulomb model.The parameter β is commonly termed as the friction angle (in case of dealing with the Drucker-Prager theory), which can be computed as the slope of the yield function in the t-p meridian plane (see Fig. 10a.By a similar computing approach for K, it is possible to derive β according to Eq. ( 24).6sin arctan( ) 3 sin It is worth mentioning that both the parameter can be defined in ABAQUS as a function of temperature and other field variables.However, they are considered as constant variables in this study because of the negligible effect of the other fields.The parameter d, called the cohesion of the material (in case of dealing with the Drucker-Prager model), is defined automatically by ABAQUS according to Eq. (25).
where, u c  is the unconfined (uniaxial) compression stress.An isotropic hardening is implemented for this model in ABAQUS by a sub-option called "Drucker-Prager Hardening" at the property module.Therefore, the tabular data from the compression post-yield stress-strain diagram can be used to define u c  (and accordingly d) as a function of the absolute plastic strain p  .The plastic flow of the linear Drucker-Prager is considered by means of a flow potential function LDP G given by Eq. (26). tan The parameter ψ, called dilation angle, specifies the direction of the plastic strain flow.As can be seen in Fig. 10a, the plastic flow is associative if ψ = β, otherwise it is non-associative.The parameters to identify the material model are reported in Tab. 6.The software automatically computes the cohesion by giving the yield stress in compression as a function of the plastic strain allowing the isotropic hardening of the yield function.In this study, the failure criterion is intended as a yield surface, i.e. when Eq. ( 21) is satisfied, plastic flow occurs.Hardening of the yield function follows until the maximum strength is reached, afterwards failure is modeled by softening.Hardening (intended as hardening followed by softening) of the yield function (see Fig. 11) is described by the equation proposed by Lubliner et al., [42] for concrete as defined in Eq. ( 27).
Eq. ( 27) expresses the flow stress  in a uniaxial test as a function of the axial plastic strain p  .The elastic properties are taken from [36].The Mohr-Coulomb friction angle is obtained through the procedure described in [43] using the properties reported in Tab. 7. To assure convexity of the yield function, the value of the flow-stress ratio is limited by a minimum value of 0.778.Since the obtained value is lower, the minimum value allowable is used instead resulting in the criterion obtained not being equivalent to the starting Mohr-Coulomb criterion but being its best approximation [30].A good approximation of the dilation angle for rock showing brittle behavior is equal to 1/4 times the friction angle [44].In Tab. 7, the parameters to define the curve of Eq. ( 27) are reported.The parameters a and b, which are obtained using Eqs.( 28) and ( 29), are determined equal to 2.68 and 665.32, respectively.It is worth mentioning that the plastic strain at Eq. ( 29) (to define the parameter b) is determined when the derivative of Eq. ( 27) is equal to zero.This plastic strain corresponds to the experimental data at which the maximum stress is reached.

Simulation of the four-point bending test
Replication of the Flexural test was obtained by means of two numerical models that consists of rollers, compressive platens, steel blocks and the specimen.Due to the high computing time required for the KCC, and by help of the symmetric nature of the test, only one-quarter of the test was modelled in LS-DYNA (see Fig. 12a).Since the ghost particle methods (which is explained below) is not yet available in ABAQUS, the numerical model replicated in this software consists of the full model (see Fig. 12b).In all the models, the conversion to SPH particles is only considered for the rock specimen.The axes of rollers are fixed in the XY plane, while the displacement-controlled compressive loading is applied by the upper platens.
The lower platens are limited to zero degree of freedom to represent the bed of the testing machine.The mechanical properties of these components are reported in Tab. 8.The hexagonal constant stress elements are used for all the solid parts.The automatic penalty-based formulation contacts are considered for all the components and the static friction coefficient is set to 0.4.In the LS-DYNA model, instead of applying single point constraints to the SPH particles, which can lead to inaccurate results and numerical instabilities, specific boundary conditions at the symmetry planes were imposed [39].The BOUNDARY_SPH_SYMMETRY_PLANE keyword creates automatically an imaginary plane which reflects the forces of a set of ghost particles to the particles in the model.Although these ghost particles have identical properties as the real ones, they do not physically exist and simply contribute to the particle approximation [21].The maximum principal strain is considered as the eroding criteria for FEM to SPH particles conversion for all the numerical simulations.The numerical results of KCC and LDP models in terms of stress distribution along length are shown in Fig. 13a and b, respectively.Fig. 13 represents the distribution of the X stress (along the length of specimen) one time step before failure.The tensile and the compressive stresses at both numerical models are present in the element below and above the neutral axis, respectively.However, the amount of maximum stresses in the two models are different.This issue is discussed in next chapter by comparison of the numerical results with the experimental ones.

COMPARISON OF THE NUMERICAL RESULTS WITH EXPERIMENTAL DATA
he load-displacement diagrams of the numerical models; developed by the KCC (both the automatic input generation mode and the calibrated mode) and the LDP material models are compared to the experimental results according to Fig. 14 and Tab. 9.There is a row in Tab. 9, labeled as "95% CI" of experimental average value, that was also reported in Tab. 1.As previous described there is a probability of 95% that the average value of a set of populations lie inside this interval.Therefore, they can be considered a reasonable description of the physical phenomena thus suitable for comparing the numerical results.The KCC material model shows the best response compared to the others, since the reported values of all three parameters, i.e.Maximum load, Maximum displacement and Flexural strength, lie within the 95% Confidence Interval of the average value.As expected, the results obtained by the automatic input mode of the KCC model are not in good agreement with the experimental data.This issue is more critical in the displacement level, and the flexural strength is very close to the minimum value reported as 95% Confidence Interval.LDP predicts the flexural strength as 16.25 MPa which is far from the reported 95% Confidence Interval (with an overestimation of about 80% respect to the experimental average value).Therefore, it is not reliable to model the mechanical behaviour of Pietra Serena in case of a significant tensile stress is present.Indeed, the Drucker-Prager criterion is well suitable to model the linear dependence of the material strength with p.The material behaviour of rocks, for negative values of p, deviates from linearity and the material failure points lie below the Drucker-Prager failure line leading to the overestimation of the flexural strength of the material.The numerical results obtained by the full calibrated mode of KCC, on the other hand, show significant agreement with the experimental data.By recalling the fact that this material model is calibrated based on the Berea sandstone (experimental data of the triaxial compression test), two main conclusions can be drawn.Apart from the authors' hypothesis of a similar mechanical response of Pietra Serena and Berea sandstone, the efficiency of the KCC model for numerical modelling of rock materials, and in particular sandstones, has been evaluated here. Maximum

CONCLUSION
he Four-Point Bending tests were performed based on the protocols of the ASTM standard.The average experimental flexural strength of Pietra Serena was measured as 8.98 MPa and Pietra Serena and Berea Sandstone show similar mechanical properties.The numerical simulations using the FEM-coupled to-SPH method were implemented in conjunction with the KCC and LDP material models.The KCC model showed a significant improvement when the material input parameters were calibrated and directly inserted into the material keyword, instead of the automatic input generator mode.While the LDP material model showed some limitation in analyses of the state of stress in which tensile behaviour is present.The specifications of the KCC calibrated material model are summarized in Tab.10.The flexural strength obtained by the numerical models developed by the LDP and the KCC automatic input mode lie significantly out of the 95% Confidence Intervals of the average value, as predicted from statistical elaboration of experimental data, while, this issue is solved by implementing the calibrated input mode of KCC.The coupling of SPH and FEM was proved to be a reliable method to simulate crack generation.This method was shown to be an effective solution to simulate several rock mechanics applications that deal with large deformations, e.g.penetration, drilling, perforation, etc.The numerical results are expected to be improved by a direct material identification based on the Pietra Serena, i.e. the triaxial compression and isotropic compression tests.

Figure 2 :
Figure 2: The extensometer set up on the flexural test configuration.

Figure 3 :
Figure 3: The broken specimens after performing the flexural test; (a) front view; (b) isometric view.

Figure 4 :
Figure 4: The load-displacement diagram of the flexural test.The shape of the load-mid span displacement, represented in Fig. 4, indicates a significant non-linear behaviour until the failure.The failure is denoted by a sudden drop of the loading level.A study by Meda et al. [32], investigated Pietra-Serena sandstone under a flexural test.They performed the tests on specimens with different shapes and considered a notch in the middle-span of the specimens to assure that the crack takes place at the center.The experimental data provided by Meda et al. are in a good agreement with the results obtained during this research work where the geometries are similar to each other; i.e. the flexural strength and the load-displacement diagrams of the Meda's work were reported as 9.1 MPa and 3.4 mm, respectively.

Figure 5 :
Figure 5: SPH particles in a 2D problem domain.
to 2000[26][27][28][29], known as the "Karagozian and Case Concrete (KCC or K&C)" model, by applying several significant enhancements to the material model 16.The Karagozian and Case concrete model, implemented in LS-DYNA, consists of three-invariant failure surface formulation, i.e. 1 2 3 plane defined by means of the above-mentioned equations is shown schematically in Fig.7.

Figure 7 :
Figure 7: Schematic representation of the KCC failure surfaces in the compressive meridian.

Figure 8 :
Figure 8: Deviatoric section proposed by Willam-Warnke model KCC model also equal to t c r r ), both the numerator and the denominator of the right-hand side of equation are divided by 2 same as the yield strength level y   at this range.A hardening plasticity response occurs after yielding and before reaching the maximum strength.The KCC material model determines the current state of stress at this range (between point 1 and point 2) according to Eq. (17).

Figure 9 :
Figure 9: a) the linear interpolation between failure surfaces, b) the damage function, c) unconfined compression stress-strain diagram.
It increases up to unity at a user-defined value m  , corresponding to point 2, where the maximum failure surface has been reached.The KCC model considers the hardening plasticity by means of this linearpiecewise function of η -λ.After point 2, where the softening takes place, η decreases to zero at end  corresponding to point tabulated damage function of the KCC model can be determined by up to 13 pairs of η -λ parameters.The variation of η and λ is summarized in Tab. 2. η λ Current failure surface position are pressure (positive in compression) and maximum tensile strength, respectively.The p d is, the effective plastic strain increment equal to(2 3)

Figure 11 :
Figure 11: The post-yield stress-strain diagram in compression.

Figure 13 :
Figure 13: Distribution of numerical stresses in the X direction (along the length); (a) KCC model; (b) LDP model.

Fig. 15
Fig.15indicates the comparison of the numerical results and the experimental test after failure.The contour plots represented in Fig.15a and brepresent the "scaled damage measure, δ" and the maximum principal plastic strain at failure, respectively.As can be seen, the SPH particles formed at the instant of failure express the crack patterns observed during the experimental test (see Fig.15c).

Figure 14 :
Figure 14: The load-displacement diagram; the numerical results compared to the experimental ones.

Figure 15 :
Figure 15: The comparison of the numerical results and the experimental test after failure.

Table 1 :
also expresses the load-mid span displacement diagram of the all specimens.The experimental results of the flexural test.

Table 2 :
The KCC damage evaluation parameters.

Table 3 :
The experimental data of Berea sandstone provided for the automatic mode.
3, means the UCS is considered equal to 62[MPa].Tab. 4 expresses the results obtained from the automatic input generation mode of MAT_072R3, after performing the initial simulation, which were also used at the full input calibration mode.

Table 4 :
The results obtained from the automatic input generation mode of MAT_072R3.

Table 5 :
The ai-parameters experimental of Berea sandstone provided for the full input mode.

Table 7 .
The experimental data of the Pietra-Serena sandstone.

Table 8 :
The mechanical properties of the rollers, platens and steel blocks.

Table 9 :
The numerical results of the flexural test, and the "95% CI" with respect to the average experimental data.

Table 10 :
The specifications of the full input mode of the KCC material model implemented for Berea sandstone.