Determination of Young’s Modulus, Fracture Energy and Tensile Strength of Refractories by Inverse Estimation of a Wedge Splitting Procedure

The wedge splitting test according to Tschegg provides a technique to characterize the fracture behaviour of ordinary ceramic refractory materials. By fitting the data from finite element simulation to the results of the wedge splitting test, Young’s modulus and parameters describing the failure behaviour under Mode I conditions can be inversely estimated through an adaptive nonlinear least-squares algorithm. The results show Young’s modulus is accurately identified as well as the tensile strength and total specific fracture energy when a trilinear strain softening law is employed. The inversely estimated parameters from three experimental curves of the same material at room temperature are very consistent as well as the values of thermal stress resistance parameter R and characteristic length l ch . The method developed enables the identification of the total specific fracture energy, tensile strength and Young’s modulus with numerically robust method in the relatively short time from a single wedge splitting procedure.


Introduction
Considerable deviations from pure linear elastic fracture mechanics have to be expected when the process zone size (the length of the process wake a plus the size of the frontal process zone, see Fig.1) is not negligibly small relative to the crack length or the specimen size [1][2][3][4]. This might be the case for materials showing a considerable size of the structural elements (e.g. large grains) relative to the size of the whole part and not behaving totally brittle. Examples are common building materials (especially concrete and mortar) and refractories. For those cases a material model allowing numerical treatment of crack propagation is desirable. For this purpose the fictitious crack model according to A. Hillerborg was already introduced 30 years ago and especially applied for concrete [1,2]. It assumes a strain softening behaviour characterized by a monotonously decreasing stress transferred between the crack faces until an ultimate crack opening of ult X is achieved (Fig.2a), whereas a normally linear stress/strain relation is applied in the uncracked region (Fig.2b). As Eq. (1) shows, the area below the strain softening curve is equal to the total specific fracture energy Gf.
The wedge splitting test according to Tschegg [5,6] is suitable for the determination of the specific fracture energy of refractories by favoring stable crack propagation in specimens with sufficiently large dimension [7]. The detailed description of the wedge splitting test and its applications are available in references [8][9][10][11][12], and here only the schematic representation of this test is shown in Fig.3a. Out of this test, a load/displacement curve can be registered (Fig.3b) and the specific fracture energy can be further obtained by calculating the area under the load/displacement curve relative to the projected area of the fracture surface. In many cases the test could be performed until the descending load approaches to zero, and thus the total specific fracture energy can be gained practically. Nevertheless, in the case of refractories with largely reduced brittleness, the wedge might eventually hit the groove of the specimen before the load approaches to zero. A premature termination of the test is necessary for safety considerations and hence only a major part of the specific fracture energy (denoted as G'f) can be calculated by integration: where ult  is the ultimate displacement and A the area of the projection of the fracture surface. As a compromise the test may be terminated in case the descending load reaches 15% of the maximum force. Besides, a nominal notch tensile strength NT is yielded which comprises both tensile stresses and flexural stresses in the ligament of the specimens and can be determined: where b and h are the width and height of the fracture surface area, respectively, and y is the vertical distance of the centre of gravity of the fracture surface from the horizontal force. Direct determination of the tensile strength ft is not possible from the wedge splitting test out of several reasons. These are the actions of the notch, the additional flexural load and the size effect which tends to reduce the apparent strength with rising specimen size due to increasing brittleness.
To define the total specific fracture energy and tensile strength, an inverse estimation procedure could be applied for the identification from the wedge splitting test results. Pioneer work was performed by A. Hillerborg [13], and P.E. Roelfstra [14] published similar studies later, in which an optimization procedure was established to fit the experimental and predicted load/displacement curves from FEM (finite element method) simulation of 3-point bending tests with a bilinear strain softening law. However, the fitting was not really satisfying especially for the area near the peaks of the curves and similar phenomena were observed in the work of J. Kim where the same inverse estimation strategy was used for the wedge splitting test [15], probably because the optimization unfortunately fell into a local minimum. In order to avoid the risk of local minima an evolutionary algorithm was employed and extended for the optimization during the inverse determination of a strain softening law [16]. However, this optimization method did not show high efficiency and an extreme long computation time was needed for one optimization.  [8][9][10][11][12].
Several analytical models or approaches for the wedge splitting test of concrete and mortar were proposed which can inversely calculate the strain softening curves rather efficiently, for instances, hinge model and back analysis [17][18][19][20][21][22][23]. A common character of them is the inverse analysis has to be performed step by step according to elastic and different stages of strain softening behaviour. Consequently the interactions among the modulus of elasticity, tensile strength and fracture energy might be ignored or inhibited. Frequently, a weighting function was also demanded in order to reach a better local and global fitting on the experimental curves [24]. Besides, the optimization algorithms used could not exhibit favorable convergence by using the common least squares approach and thus a specific error norm or peculiar convergence tolerance has to be adopted [20,25]. In the refractory field, trials also have been done to inversely determine Young's modulus E as well as the tensile strength and fracture energy from the wedge splitting test [26]. The FEM simulations of the wedge splitting test were performed in terms of a strain softening model according to Hordijk [27], and the displacement of the descending branch at 50% of Fmax was chosen as a significant feature accounting for the shape of the load/displacement curve. Three nonlinear formulas were thus established to numerically describe the relations of fracture mechanical properties of refractories. By a minimization technique, the tensile strength, specific fracture energy and Young's modulus have been inversely estimated.
In the present paper a more general approach is shown to inversely identify the tensile strength, the total specific fracture energy and Young's modulus of refractories from wedge splitting test results. By this approach, a pre-analysis of Young's modulus from the experimental curve is not necessary, and the definition of the stress/strain softening type can be easily done in the commercial code ABAQUS [28]. Moreover, the conventional leastsquares norm is employed for the minimization of the difference between experimental and predicted data. The load/displacement curves of a magnesia-chromite material determined by the wedge splitting test at room temperature are used for a sensitivity investigation of number and values of fracture mechanical properties. The final inverse estimation results of fracture mechanical properties are provided and compared to experimental results.

Mathematical procedure
An inverse estimation procedure usually can be realized by minimizing the objective function through searching a better parameter vector. Frequently, the objective function ) (x f is given as the sum-of-squares as follows: is the residual vector, n is the number of its components, x is the parameter vector, and p is the number of parameters.
Several algorithms such as Gauss-Newton and Levenberg-Marquardt methods were developed to solve the least-squares problems and well applied for inverse estimation [29][30][31]. NL2SOL, an adaptive nonlinear least-squares algorithm, is more promising in dealing with large residue problems i.e. least-squares problems where the residuals do not tend to zero with increasing number of iterations [32,33]. In this algorithm, the increment of x for iteration is determined customarily by Here g is the gradient of where J is the Jacobian matrix of ) (x W . H is defined according to following two equations: Eq. (9) represents the Gauss-Newton model (abbreviated as G) and Eq. (10) the augmented model (denoted as S) which additionally includes the Hessian matrix The choice of the model is decided by evaluating the decrease of the actual function result relative to its predicted value at each iteration and this strategy can promote the fast convergence even the starting guesses are poor.

Models and inverse procedure
A two-dimensional and symmetrical model comprising one half of the specimen and the wedge as well as one transmission piece was built according to the test configuration and specimen geometry (Fig.4a). Strain softening behaviour was assigned to quadrilateral cohesive elements with linear interpolation scheme arranged along the ligament area which was 1.5 mm wide (shaded area in Fig.4b) and had a projected fracture surface of 65X63 mm 2 . Linear elastic behaviour was assumed to the bulk of the material and the load transmission pieces with quadrilateral elements interpolated by linear scheme. The modeling of the wedge splitting test process was performed via the commercial code ABAQUS. Two executable files were written in Python language in order to extract data from ABAQUS .odb file as well as discretize the experimental data. The open source code DAKOTA in which the algorithm NL2SOL has been implemented was used to integrate all the components of the minimization procedure [34]. The numerical process is shown in Fig.4b.

Sensitivity investigation
A magnesia-chromite material (56.5wt% MgO, 25.5wt% Cr2O3) was chosen for the wedge splitting test at room temperature. Three specimens were prepared according to the above mentioned geometry and load/displacement curves were recorded. Two hundred equally spaced data points were extracted from each experimental curve, which were believed to contain sufficient information of the original load/displacement curve. For the sensitivity investigation, only one curve was adopted and six cases of parameters were designed for the evaluation as seen in Table 1. In the case of an odd number in the first column Young's modulus was one of the parameters identified. On the other hand, in the cases with even number Young's modulus was measured by Resonance Frequency and Damping Analyzer (RFDA, IMCE in Belgium) at room temperature, which is 81 GPa. The corresponding type of strain softening laws also can be seen in the last column. Table1 Testing designs for sensitivity investigation.   Fig.5 shows an example of the convergence and alternative numerical models during the minimization. Herein the relative residual was defined as the difference between the previous and the current residual results divided by the previous residual value [33]. G-S in the figure means that during the iteration the Gauss-Newton model was tried first and a switch was then made to the augmented model. S-G and G-S-G have analogous meanings. Clearly the residual decreased dramatically during the first 3 iterations and the relative residuals were close to 1. Afterwards, the residual values changed slightly and consequently the relative residual decreased by a power of ten. The Gauss-Newton model was chosen primarily until the relative residual was 0.0015, and then the augmented model was chosen to enhance the convergence where the relative residual was increased to 0.004. Afterwards, the switch between G and S models was more frequently among the 7 th -10 th iterations. The minimization procedure was automatically terminated when the relative residual was less than 10 -6 , and a further decrease of the residual was not possible by this method.
The final residuals and number of function evaluations were used to compare the convergence capacity of the inverse estimation. The number of function evaluations is the number of calls to calculate the objective function ) (x f , which includes those used for the evaluation of Eq. (2). As can be seen from Fig.6 if more than four parameters are evaluated the residual is only the seventh part of that in the case with two parameters. In contrast, the number of function evaluations increases only by a factor of two or three. Usually, a typical computation time was around 0.5 h for each case.
The influences of Young's modulus on the identified results for the total specific fracture energy Gf and the tensile strength ft were analyzed as well. Figures 7-8 show normalized values for Gf, ft and E. In the case of Gf and ft they are obtained by dividing the actual value by the maximum occurring in the test series, and for E the measured Young's modulus by RFDA is the reference value. The values attributive to the dots in Figs.7 and 8 indicate the number of evaluated parameters. It is shown that Gf did not receive significant impact from the variation of Young's modulus. Nevertheless, the normalized ft varied from 0.83 to 1. Because the strategy is to minimize the difference of two curves and Gf is the specific area under the curve, Gf did not vary too much although Young's modulus changed. Only cases containing 4, 6 and 7 evaluated parameters exhibited very similar values. Frequently, thermal stress resistance parameter R and characteristic length lch (proportional to the R'''' parameter and inversely proportional to a brittleness number) were used to assess the resistance of refractories to crack initiation and propagation due to thermal shock [12,35]. These two values can be calculated utilizing the results from the wedge splitting test, Poisson's ratio (v) and thermal expansion coefficient (α). Fig.9 depicts the influence of the number of parameters on R and lch.
The maximum values of R and lch were used as denominators for the normalization. Similar to Fig.8, the application of 4, 6 and 7 evaluated parameters resulted in nearly the same values of R and lch. Although for 6 and 7 evaluated parameters lower residuals occurred, confident results for the characterization of thermal stress resistance can be obtained by a bilinear strain softening law when the Young's modulus was determined by RFDA. In the case that Young's modulus is not available, a trilinear strain softening law can be used to determine Young's modulus, the tensile strength and the total fracture energy accurately.

Results
Four parameters as shown in Table 1 were used for the inverse estimation of 3 experimental curves. Curves 2 and 3 were shifted several millimeters in order to be better distinguished. Fig.10 shows the good agreement and the residuals between the experimental and simulated curves. Here an example was provided to show the identified parameters of bilinear strain softening law from Curve 1: ft=5.50 MPa, σ1=2.93 MPa, X1=22.2 m, Xult= 80.0 m. Table 2 lists the calculated and inversely estimated parameters from the 3 experimental curves and their standard deviations. It is shown that the mean value of ft equaled 5.5 MPa while that of NT was 8.9 MPa. A difference between NT and ft is expected due to not only the flexural loading component in NT but also the size effect. It is commonly recognized that the nominal tensile strength increases with decreasing structural size in the quasi-brittle materials especially in the presence of notches or cracks [36][37][38]. Several general laws comprising the structural sizes and characteristic length (or brittleness number) were developed to describe    Residual /N 2 Normalized lch [1] the size effect on the nominal tensile strength in the last decades [39][40][41][42]. Especially, a formula of the ratio NT/ ft dependent of the brittleness number was given in the case of wedge splitting test according to Tschegg, and this ratio decreased with increasing brittleness number [26]. In the case of concrete or mortar with large brittleness, the difference between the nominal notch tensile strength and the tensile strength is expected to decrease.
As seen in table 2, the mean value of G'f was equal to151 Nm -1 and that of Gf 177 Nm -1 which was 17% higher than the former. This difference between G'f and Gf is reasonable because G'f was calculated up to a fixed level of the residual load in laboratory as mentioned early. Also, R and lch were calculated from the above determined parameters respectively and listed in Table  2. The mean value of R derived from the measurement was 8.3 K and that from the inverse estimation only 5.2 K. Evidently difference between the values of lch derived from the inverse estimation and the measurement can be observed. Normally the characteristic length is based on a tensile strength and therefore the value calculated by the results of the proposed method is preferred. In other words, the real lch value was underestimated by 3 times if lch was calculated from the experimentally determined parameters. This is largely influenced by the fact that the nominal notch tensile strength is considerably larger than the identified tensile strength. Additionally, the inverse estimation results showed comparatively small deviations for the three curves although their residuals were apparently different (Fig.10). A strong robustness of this inverse estimation procedure can be expected for the wedge splitting test.

Conclusions
A computational procedure has been successfully applied for the wedge splitting test in order to inversely estimate not only strain softening law under Mode I failure but also Young's modulus of ordinary refractory ceramic materials. The inverse estimation results from three experimental curves of one material show the good agreement. The tensile strength largely differs from the nominal notch tensile strength, which means that this inverse evaluation procedure is very necessary to do if the strength should be assessed. The difference in the specific fracture energy is by far smaller. In the case of table 2, G'f is 85% of Gf, but NT is 162% of ft. It is worthy to comment that this general inverse estimation procedure is easily handled and may be implemented into other research fields where the inverse estimation is needed.