Improvement of Mixed-Mode I/II Fracture Toughness Modeling Prediction Performance by Using a Multi-Fidelity Surrogate Model Based on Fracture Criteria

Today, artificial intelligence plays a huge role in the mechanical engineering field for solving many complex problems and the problem with fracture mechanics is one of them. In fracture mechanics, artificial intelligence is used to predict crack behavior under various conditions such as mixed-mode loading. Many parameters are used for explaining the crack behavior under various conditions, but those parameters are obtained from destructive testing, in which usually, only one data point is obtained from each test. An artificial problem method requires a large amount of data to train the model to be able to learn crack behavior, which is a disadvantage of applying this method to fracture mechanics. To eliminate the disadvantage of the large amount of experiment data required for modeling, in this study, the small data obtained from the experiment along with data obtained from fracture criteria that were used for elementary prediction of mixed mode fracture toughness were used to create an artificial intelligence model. Data from the experiment was combined with fracture criteria data using the multi-fidelity surrogate model that is described in this study. The mixed mode I/II fracture toughness of the PMMA material was tested in order to primarily propose the data combination technique. After the modeling process, the prediction results indicated that the performance of a model in which the actual test data was combined with the data from the fracture criteria (multi-fidelity surrogate model) was more predictively effective compared to only actual data-based modeling.


Introduction
One of the crucial factors to consider when designing engineered components is fracture toughness. This parameter of the fracture mechanics field indicates a material's capacity to withstand the crack expansion of discontinuous or cracked materials in the presence of external loads. According to the direction of the external load acting on the crack surface, fracture toughness can be divided into three modes [1]. For the majority of engineering parts, the load pattern or direction is more likely to occur in a mixed direction than in just one direction. According to the nature of the load direction on that part, it must be studied in a mixed-mode fracture toughness; either it is studied in mixed-mode I/II [2,3], mixed-mode I/III [4,5], or mixed-mode II/III [6,7], with mixed-mode I/II loading being popularly studied in various research. The fracture toughness values under mixed-mode loading can be calculated under several methods such as actual testing, numerical method, etc. In addition to such methods, in recent years, the artificial intelligence method was used to calculate or predict fracture toughness values.
Artificial Intelligence is one of the most effective solutions to complex problems based on human learning and decision-making. The majority of artificial intelligence models that predict fracture toughness use supervised learning which trains their models using data from actual fracture toughness tests in order to produce accurate predictions [8][9][10]. For this reason, the dataset used to train the model must be sufficiently large enough so that the model can learn the behavior of the material. The relatively data-intensive use of modeling artificial intelligence is considered a huge disadvantage in fracture toughness prediction. The fracture toughness test is destructive testing in which the test specimens are destroyed during testing and cannot be used for other testing anymore. According to the testing characteristics of fracture toughness, tests that are performed in several types of engineering materials such as polymers, metal, composites, or building materials usually obtain a single point of data from each test. As a result, gathering sufficient data to model artificial intelligence necessitates extensive actual testing, which appears to go against the idea that predictive modeling is used to minimize testing. In addition to the fracture toughness data from the actual test, when considering the method for determining the fracture toughness, a method by which the elementary fracture toughness can be calculated is known as fracture criteria, whereby a lower cost of data acquisition is compared to the actual test data.
The fracture criteria is a numerical model based on the relationship between stress, strain, or energy occurring around the crack tip of material. The fracture criteria was designed for pure model loading, but it can be extended to mixed-mode loading. Ordinarily, the fracture criteria is designed to indicate whether material is fractured or not by giving a specific value, but it can also be used as an elementary prediction of fracture toughness. Many fracture criteria are created based on the relationship around the crack tip. as Also of interest is the generalized maximum tangential stress criteria (GMTS) [11,12], which focuses on the tangential stresses at the crack tip, and the average strain energy density criteria (ASED) [13,14], which focuses on the strain energy density around the crack tip, and the generalized maximum energy released rate criteria (GMERR) [15], which concerns the strain energy release rate around the crack tip etc. The fracture criteria serves as an excellent indication of the fracture of the material, but sometimes the specific value of fracture toughness obtained from the fracture criteria is quite different from the actual fracture toughness, as seen in the previous work of Poapongsakorn et al. [5], which describes the time-dependence of epoxy resin that was found in the elementary fracture toughness from criteria clearly different from the actual fracture toughness value.
In order to solve the relatively data-intensive use of modeling artificial intelligence, as in the case of models where the fracture toughness of materials is applied and where the acquisition of data is limited due to high cost or limited material, or when it is difficult to prepare specimens such as composites materials that take a lot of time to be molded or metal materials that take a long time in the heat treatment process, etc., this research aims to create an artificial intelligence model using small data from actual fracture toughness tests. The data from the actual fracture toughness test are combined with the elementary prediction data obtained from the widely used fracture criteria, including the generalized maximum tangential stress criteria (GMTS), the average strain energy density criteria (ASED), and the generalized maximum energy released rate criteria (GMERR) in the model training stage to obtain enough data to train the model to be able to make accurate predictions. Poly (methyl methacrylate) or PMMA, which is frequently used in lab-scale experiments, was chosen to perform fracture toughness testing under mixed-mode I/II loading in order to concentrate on presenting an artificial intelligence modeling technique from a small actual test dataset by combining the data with fracture criteria. The incline crack specimen with three-point bend testing configurations was selected as described in the first section of this study. The widely used fracture criteria for elementary fracture toughness prediction are described in the second section. Combining actual testing data with elementary prediction data from fracture criteria for artificial intelligence modeling relies on the multi-fidelity surrogate model integration principle [16][17][18] commonly used for multiple data combinations in aircraft design and the Kriging algorithm [19][20][21], which is designed for modeling artificial intelligence from small datasets that are described in the third section. The performance of the artificial intelligence model generated from the data integration is compared with the performance of the model generated from the only actual test data shown in the last section of this study.

Mixed-Mode I/II Fracture Toughness
This study aims to present a technique for improving artificial intelligence model prediction accuracy using a limited testing dataset. A poly (methyl methacrylate) or PMMA sheet, which is widely used in lab-scale testing, was employed to test the fracture toughness under mixed-mode I/II in order to concentrate largely on the presentation of the approach. Table 1 shows the mechanical properties of the PMMA sheet utilized in this study, for which testing was conducted in accordance with ASTM D638 and D732 standards. Numerous testing methods and specimens, such as incline edge crack asymmetric bending [3], compact tension shear [6], asymmetric four-point bending [22], semi-circular bending [23], Brazilian disc specimen [24], etc., are used to study mixed-mode I/II fracture toughness. In this study, the mixed-mode I/II of PMMA was studied using the inclined crack bending specimen based on previous research by Mingdong Wei [25] and M.R.M. Aliha [26]. The inclination crack bending specimen (ICB) was modified to be more suitable for our three-point bending testing apparatus by extending its width and length. The dimension of the ICB specimen is shown in Figure 1 (the crack length, inclined angle, span length (S), and thickness of specimen (t) will be explained in the next section).
Materials 2022, 15, x FOR PEER REVIEW 3 of 23 toughness prediction are described in the second section. Combining actual testing data with elementary prediction data from fracture criteria for artificial intelligence modeling relies on the multi-fidelity surrogate model integration principle [16][17][18] commonly used for multiple data combinations in aircraft design and the Kriging algorithm [19][20][21], which is designed for modeling artificial intelligence from small datasets that are described in the third section. The performance of the artificial intelligence model generated from the data integration is compared with the performance of the model generated from the only actual test data shown in the last section of this study.

Mixed-Mode I/II Fracture Toughness
This study aims to present a technique for improving artificial intelligence model prediction accuracy using a limited testing dataset. A poly (methyl methacrylate) or PMMA sheet, which is widely used in lab-scale testing, was employed to test the fracture toughness under mixed-mode I/II in order to concentrate largely on the presentation of the approach. Table 1 shows the mechanical properties of the PMMA sheet utilized in this study, for which testing was conducted in accordance with ASTM D638 and D732 standards. Numerous testing methods and specimens, such as incline edge crack asymmetric bending [3], compact tension shear [6], asymmetric four-point bending [22], semi-circular bending [23], Brazilian disc specimen [24], etc., are used to study mixedmode I/II fracture toughness. In this study, the mixed-mode I/II of PMMA was studied using the inclined crack bending specimen based on previous research by Mingdong Wei [25] and M.R.M. Aliha [26]. The inclination crack bending specimen (ICB) was modified to be more suitable for our three-point bending testing apparatus by extending its width and length. The dimension of the ICB specimen is shown in Figure 1 (the crack length, inclined angle, span length ( S ), and thickness of specimen ( t ) will be explained in the next section).  Stress intensity factors ( K ) were used to present the mixed-mode I/II fracture toughness values of PMMA in the form of the ICB specimen shape, which were divided into mode I ( I K ) and mode II ( II K ) stress intensity factors, which we will refer to as "fracture toughness" to avoid confusion. The I K and II K of the ICB specimen can be calculated using the Equations (1) and (2), respectively. For ICB specimens, the crack's inclination at an angle to the line of load applied to the specimen causes a mixed load  Stress intensity factors (K) were used to present the mixed-mode I/II fracture toughness values of PMMA in the form of the ICB specimen shape, which were divided into mode I (K I ) and mode II (K I I ) stress intensity factors, which we will refer to as "fracture toughness" to avoid confusion. The K I and K I I of the ICB specimen can be calculated using the Equations (1) and (2), respectively. For ICB specimens, the crack's inclination at an angle to the line of load applied to the specimen causes a mixed load between the open crack loading (mode I) and the shear crack loading (mode II) to be created. The ratio between mode I and mode II loading, also called "mode mixity parameter (M e )", can be calculated using Equation (3) where P is maximum load applied to the specimen, W is the specimen width, t is the specimen thickness, a 0,ICB is the initial crack length, S is the span length that measures the distance from the center of the specimen to the bottom of the three-point bending pin in horizontal axis, β is the inclined crack angle, Y I and Y I I are the normalized or dimensionless stress intensity factors for mode I and mode II, respectively, which depend on the specimen's geometry and loading configurations. The corresponding values of Y I and Y I I of the ICB specimen at different specimen geometries and loading configurations can be computed using the finite element analysis model (FEA). The FEA model of the ICB specimen in ABQUS 6.13 software is shown in Figure 2. The FEA modeling used material properties of PMMA (Table 1) and the load applied to the model was fixed at 1 kN. The FEA model of all specimen geometries and loading configurations was created by twenty nodes quadratic C3D20 elements (Hex-dominated) with a convergence at 1 mm, the approximate global element size. Mesh refinement was generated around the crack tip for singularity concerns of the stress-strain field around the crack tip ( Figure 2 right-hand side) by the fifteen nodes quadratic C3D15 element (Wedge) with a 0.15 mm element size. For elements that are defined according to the above, the model had approximately 21,536 elements. According to actual testing, the moving condition of the FEA model was set, and the support pin on the bottom left side is set to be motionless or fixed to both move and rotate. The support pin on the bottom right side was configured as a roller support that can be moved along the x-axis and rotated in the direction normal to the x-y plane.
where P is maximum load applied to the specimen, W is the specimen width, t is the specimen thickness, 0,ICB a is the initial crack length, S is the span length that measures the distance from the center of the specimen to the bottom of the three-point bending pin in horizontal axis, β is the inclined crack angle, I Y and II Y are the normalized or dimensionless stress intensity factors for mode I and mode II, respectively, which depend on the specimen's geometry and loading configurations. The corresponding values of I Y and II Y of the ICB specimen at different specimen geometries and loading configurations can be computed using the finite element analysis model (FEA). The FEA model of the ICB specimen in ABQUS 6.13 software is shown in Figure 2. The FEA modeling used material properties of PMMA (Table 1) and the load applied to the model was fixed at 1 kN. The FEA model of all specimen geometries and loading configurations was created by twenty nodes quadratic C3D20 elements (Hex-dominated) with a convergence at 1 mm, the approximate global element size. Mesh refinement was generated around the crack tip for singularity concerns of the stress-strain field around the crack tip ( Figure 2 right-hand side) by the fifteen nodes quadratic C3D15 element (Wedge) with a 0.15 mm element size. For elements that are defined according to the above, the model had approximately 21,536 elements. According to actual testing, the moving condition of the FEA model was set, and the support pin on the bottom left side is set to be motionless or fixed to both move and rotate. The support pin on the bottom right side was configured as a roller support that can be moved along the x-axis and rotated in the direction normal to the x-y plane. According to preliminary experiments, the length of the FEA model was set at 100 mm, the width at 30 mm, the ratio between initial crack length and specimen width set at 0.6, the ratio between span length and specimen width set at 0.6, and the incline crack angle ( β ) was increased from 0° to 50° (increment 5°). According to Equation (4), I Y and II Y can be converted from I K and II K extracted from the finite element analysis using According to preliminary experiments, the length of the FEA model was set at 100 mm, the width at 30 mm, the ratio between initial crack length and specimen width set at 0.6, the ratio between span length and specimen width set at 0.6, and the incline crack angle (β) was increased from 0 • to 50 • (increment 5 • ). According to Equation (4), Y I and Y I I can be converted from K I and K I I extracted from the finite element analysis using the J-integral method [26,27]. The dimensionless stress intensity factor Y I and Y I I at various incline crack angles are shown in Figure 3. The results showed that Y I decreased as the incline crack angle increased, whereas Y I I increased until reaching the maximum value at an incline crack angle around 20 • and then decreased. According to the FEA result, the incline crack angle was chosen to be 0 • for pure mode I (M e = 1), 20 • for mixed-mode I/II (M e = 0.5), and 31 • for pure mode II (M e = 0). Fracture toughness testing was carried out on a Lloyd universal testing machine LD series (100 kN) at a loading rate of 10 mm/min. The ICB specimen was prepared using a laser cutting process on a PMMA sheet. A brand-new razor blade was used to cut the crack tip of the ICB specimens to ensure that it contained the theoretically sharp crack. In Equations (1) and (2), the mixed mode fracture toughness shown values were changed according to the function of stress (P/Wt) applied to the specimen. The stress function applied to the specimen is found to be related to the force exerted on the specimen (P) and the cross-sectional area of the stressed area (Wt). In this study, the thickness of the ICB specimen which influenced the cross-sectional area of the stressed area was included to factor in the design of the experiment (DOE). The general full factorial design was generated, with at least three replicates in each condition and the average fracture toughness of PMMA from testing used for the present predictions method is shown in Table 2.
A brand-new razor blade was used to cut the crack tip of the ICB specimens to ensure tha it contained the theoretically sharp crack. In Equations (1) and (2), the mixed mode fracture toughness shown values were changed according to the function of stress ( P Wt ) applied to the specimen. The stress function applied to the specimen is found to be related to the force exerted on the specimen ( P ) and the cross-sectional area of the stressed area (Wt ). In this study, the thickness of the ICB specimen which influenced the cross-sectional area of the stressed area was included to factor in the design of the experiment (DOE). The general full factorial design was generated, with at least three replicates in each condition and the average fracture toughness of PMMA from testing used for the present predictions method is shown in Table 2.
(4) Figure 3. Mode I and II dimensionless stress intensity factor at various incline crack angles.

Mixed-Mode I/II Fracture Criteria
For the branch of fracture mechanics that studies the behavior of cracks under various loading types, such as crack propagation, crack severity parameters, and so on, the occurrence behavior of cracks in a material can be explained by the relationship of stress, strain, or energy exerted on the crack of the material. As a result, many researchers attempted to predict crack behavior, such as whether fractures occur on materials or not when subjected to various loading conditions, using a variety of numerical equations based on a relationship of stress, strain, energy, or material behavior known as "fracture criteria". A fracture criterion defines how fracture occurs on material by giving a specific value which could be used as an elementary fracture toughness predictive tool. Most fracture criteria are typically designed for pure mode loading conditions, but they have also been developed and applied to mixed-mode loading. In this study, the generalized maximum tangential stress (GMTS), the average strain energy density (ASED), and the generalized maximum energy release rate criteria (GMERR) were used for the elementary prediction of the mixed-mode I/II fracture toughness of PMMA. The brief fracture criteria used in this study is described in this section.

Generalized Maximum Tangential Stress Criteria (GMTS)
The generalized maximum tangential stress (GMTS) is a fracture criterion that was modified from maximum tangential stress (MTS) to improve the predictive performance of the fracture criteria. The GMTS criteria is a criterion based on the stress intensity factor (K) with added the T-stress term into the original MTS criteria. For the GMTS fracture criteria aimed at studying the tangential stresses at the crack tip (σ θθ ), these stresses can be written according to the Equation (5) [25] where r is the distance to the crack tip in polar coordinates, θ 0 is initial crack growth angle measured from the initial crack direction, T is a non-singular term or knows as T-Stress, and H.O.T. is the high order term in the stress solution. According to the GMTS fracture criteria, a material is fractured when the tangential stress at the crack tip exceeds the stress that the material can withstand, with the fracture growing in the direction of the maximum tangential stress at the crack tip. The direction of crack growth or initial crack growth angle (θ 0 ) can be described in terms of dimensionless stress intensity factors that are the basis of the fracture criteria in Equation (6) (7) and (8) where r C is the critical distance (assume specimen failure under plane strain conditions), and, thus, r C can be estimated from Equation (9) [28]) and T * is the normalized form of T-stress, which is influenced by the same factors that influence dimensionless stress intensity factors. The T-stress can be calculated using the finite element analysis method and converted to a normalized form according to Equation (10). Figure 4 shows the T * at various incline crack angles. T * is a negative value at slightly inclined crack angles (less than 15 • ) and then increases with incline crack angles.

Average Strain Energy Density Criteria (ASED)
The average strain energy density criteria (ASED) was first presented by Lazzarin e al. [29]. The ASED criteria focuses on strain energy density around the crack tip. In this criterion, the material will be fractured when an average strain energy density over a control volume reaches the critical strain energy ( C W ), which is dependent on materia properties and notch geometry. In the ICB, where the crack tip is sharp, the contro volume is a circle of radius ( C R ) with a center at the crack tip. The result of the stress based expression can be summarized together with the fracture of the aforementioned material and can be rewritten in terms of mode I and mode II local strain density following Equation (11)

Average Strain Energy Density Criteria (ASED)
The average strain energy density criteria (ASED) was first presented by Lazzarin et al. [29]. The ASED criteria focuses on strain energy density around the crack tip. In this criterion, the material will be fractured when an average strain energy density over a control volume reaches the critical strain energy (W C ), which is dependent on material properties and notch geometry. In the ICB, where the crack tip is sharp, the control volume is a circle of radius (R C ) with a center at the crack tip. The result of the stress-based expression can be summarized together with the fracture of the aforementioned material and can be rewritten in terms of mode I and mode II local strain density following Equation (11) [2] where W 1C and W 2C are the critical strain energy densities obtained from pure mode I and pure mode II loading, respectively, and can be calculated following Equations (12) and (13). W 1 and W 2 are the strain energy density generated from the applied external load. The values of W 1 and W 2 can be rewritten in terms of stress intensity factor, which is considered an elementary prediction of fracture toughness of the ASED criteria (K I,ASED and K I I,ASED ) following Equations (14) and (15) where e 1 and e 2 are given by Lazzarin et al. [29] as 0.1186 and 0.3332, respectively, K I,ASED and K I I,ASED are specific values of mixed-mode I/II fracture toughness obtained from the ASED criteria, λ 1 and λ 2 are mode I and II Williams' eigenvalues, which are equal to 0.5 [30], and R 1C and R 2C are mode I and mode II radii under plane strain conditions, which can be calculated following Equations (16) and (17).
where K I IC is mode II critical fracture toughness.

Generalized Maximum Energy Release Rate Criteria (GMERR)
The generalized maximum energy released rate criteria was proposed by Hou, Cheng, et al. [15]. The GMERR was developed by adding terms of non-singular stress based on the maximum energy released rate criteria (MERR), which predicts the initiation of the mixed mode crack by calculating the energy release rate around the crack tip. In the GMERR criteria, materials will be fractured when values of strain energy release rate reach the maximum strain energy release rate of the material. The initial crack growth angle according to the GMERR can be calculated as follows where the variables are the same according to the Equation (7) and coefficients (C) can be calculated as follows The elementary prediction fracture toughness of the GMERR criteria separated by mode I and II loading (K I,GMERR and K I I,GMERR ) can be expressed as where coefficients (B) can be calculated as follows The fracture criteria mainly require the mechanical properties or critical mode I fracture toughness of materials for elementary prediction fracture toughness values that cause materials to fracture. In this study, the parameters used for the fracture criteria equations are shown in Tables 1 and 2, and Figures 3 and 4. The specific mixed-mode I/II fracture toughness from the criteria is calculated and separated based on the aforementioned influence of ICB specimen thickness on the cross-sectional area of the stressed area shown in Figure 5. The results showed good prediction performance at pure mode I and mode II fracture toughness for both thicknesses, but at mixed-mode I/II conditions, the predictions were quite different from the experimental results, indicating a decrease in the fracture criteria's predictive efficiency. To demonstrate the fracture criteria's predictive performance more clearly, the fracture criteria's performance is measured using the MAPE performance metric, which is commonly used in regression problems according to Equation (23 [31], it was found that the predictive performance of the all fracture criteria at mode I (K I ) was relatively good (10% ≤ MAPE ≤ 20%). However, when mode II (K I I ) is considered, prediction performance is a rational prediction, which means that only data trends (20% ≤ MAPE ≤ 50%) can be predicted. Although elementary predictions with fracture criteria show some deviation from actual fracture toughness, the specific fracture toughness from these fracture criteria is still used as a preliminary criterion that determines the limit stress that materials can withstand before fractures occur, which is very useful in designing engineering parts where fracture toughness is a consideration.
where A i is the ith experiment data,Â i is the ith prediction data, and n is the number of observations.
Materials 2022, 15, x FOR PEER REVIEW 10 of 23 fracture toughness from the criteria is calculated and separated based on the aforementioned influence of ICB specimen thickness on the cross-sectional area of the stressed area shown in Figure 5. The results showed good prediction performance at pure mode I and mode II fracture toughness for both thicknesses, but at mixed-mode I/II conditions, the predictions were quite different from the experimental results, indicating a decrease in the fracture criteria's predictive efficiency. To demonstrate the fracture criteria's predictive performance more clearly, the fracture criteria's performance is measured using the MAPE performance metric, which is commonly used in regression problems according to Equation (23 prediction performance is a rational prediction, which means that only data trends (20% ≤ MAPE ≤ 50%) can be predicted. Although elementary predictions with fracture criteria show some deviation from actual fracture toughness, the specific fracture toughness from these fracture criteria is still used as a preliminary criterion that determines the limit stress that materials can withstand before fractures occur, which is very useful in designing engineering parts where fracture toughness is a consideration.
(a) Thickness 5 mm (b) Thickness 10 mm where i A is the ith experiment data,  i A is the ith prediction data, and n is the number of observations.

Artificial Intelligence Method
In recent years, many researchers have applied artificial intelligence to many complex engineering problems in many disciplines. For fracture mechanics, artificial

Artificial Intelligence Method
In recent years, many researchers have applied artificial intelligence to many complex engineering problems in many disciplines. For fracture mechanics, artificial intelligence has been applied to describe the behavior of cracks, such as in the prediction of the pure mode loading fracture parameters [8,32,33] or mixed-mode loading fracture parameters [22,34]. Artificial intelligence methods have significant advantages in terms of making accurate predictions on complex problems, but there is a rather problematic disadvantage about the learning process of artificial intelligence models. The artificial intelligence learning process simulates how the human brain responds to external stimuli and learns to respond to that stimulus in the future, and learning in this manner requires a large amount of data. The fracture toughness test is a destructive test in which one data point is obtained, which means that applying artificial intelligence predictions requires a large number of tests, which appears to be counterproductive to the predictive equation's goal of reducing the number of tests as much as possible.
In this study, the artificial intelligence model's prediction performance is improved when data from the fracture toughness test is limited. The actual fracture toughness test data is then combined with the infinite elementary predictive data obtained from the fracture criteria. The concept of the data combination technique called the "multi-fidelity surrogate model" is described in the next section.

Concept of Multi-Fidelity Surrogate Model
The concept of the combination technique is explained in Figure 6 using continuous functions with nonlinear relationships. (It should be noted that continuous functions with nonlinear relationships are unrelated to this study, as they are simply an illustration to help explain the concept of combination). The high fidelity function was used to represent the fracture toughness data obtained from actual testing. In addition, it has high accuracy but is expensive or time-consuming to obtain. The low fidelity function was used to represent the fracture toughness obtained from the elementary prediction of the fracture criteria where data acquisition is quick and inexpensive, even though accuracy may not be very high. The difference between the prediction results and the actual value of the problem is clearly shown in Figure 6 (left-hand side), but if the relationship between these two values is considered, the predicted value can be equal to the actual value of the problem according to Equation (24) when the difference is known. The actual value and prediction value difference parameters may be linear or nonlinear relationships, depending on the characteristics of the problem. This correlation occurs in such a way that the data is interchangeable and is performed on a model called the "surrogate model'. The interchangeability of the data means that rather than using an excessive amount of high fidelity data to build a model, low fidelity data can be used in conjunction with the difference parameters to build a predictive model of high fidelity data, as shown in Figure 6 (right-hand side). This model is referred to as a "multi-fidelity surrogate model".
where F(X) HF,i is the output of the high-fidelity function at index ith, F(X) LF,i is the output of the low-fidelity function at index ith, and D i is difference between the output of the high and low fidelity functions.
In the multi-fidelity surrogate model, data substitution is performed on the basis of a Gaussian process regression model known as the Kriging model. Although it uses a model to make the prediction, the Kriging model's prediction does not directly predict the difference parameters but predicts that parameter along with the value of the high fidelity data or actual data. Prior to describing how the original Kriging model and the multi-fidelity surrogate model were implemented in the fracture toughness prediction, the preparation of the data and the performance evaluation of the artificial intelligence prediction model are described in the following section.
depending on the characteristics of the problem. This correlation occurs in such a way that the data is interchangeable and is performed on a model called the "surrogate model'. The interchangeability of the data means that rather than using an excessive amount of high fidelity data to build a model, low fidelity data can be used in conjunction with the difference parameters to build a predictive model of high fidelity data, as shown in Figure  6 (right-hand side). This model is referred to as a "multi-fidelity surrogate model".  Figure 6. Concept of surrogate model for multi-fidelity function.

Data Preparation and Model Performance Evaluation
Data preparation is one of the processes that affects the predictive performance of the model. In this study, the artificial intelligence model used specimen thickness (t) and the mode parameter (M e ) as inputs and average mixed-mode I/II fracture toughness (K I and K I I ) as outputs. Six data points were used in the artificial intelligence modeling process. To reduce the different lengths of input factors in the artificial intelligence model (thickness and mixed mode parameter), input data was normalized to range 0 to 1 using Equation (25). Because of the small dataset of this study, the train and test dataset in the modeling process was selected by the holdout cross-validation method. All artificial intelligence models were created using MATLAB programing. For fracture toughness in which the data is continuous, the regression artificial intelligence model was used, in which many models' performance metrics are evaluated. The performance of the prediction model was assessed in this study using the common regression performance metrics of mean absolute percent error (MAPE), root mean square error (RMSE), and coefficient of determination (R 2 ), which can be calculated in the following Equations (23), (26) and (27). (25) where X i is the input data of the artificial intelligence model at index ith, and X max and X min are the maximum and minimum data in the input dataset.
where A is the average of experimental data and other variables according to Equation (23).

Brief of Kriging Model
The Kriging model or original Kriging model represents the single fidelity artificial intelligence model that used only experiment fracture toughness in the modeling process. In this section, the original Kriging abbreviation, also referred to as "Or-K", was discussed. The Or-K modeling process is shown in Figure 7. The prediction model of the Or-K model can be expressed as follows [16,20].
The Kriging model or original Kriging model represents the single fidelity artificia intelligence model that used only experiment fracture toughness in the modeling process In this section, the original Kriging abbreviation, also referred to as "Or-K", was discussed. The Or-K modeling process is shown in Figure 7. The prediction model of the Or-K model can be expressed as follows [16,20].
The Or-K model prediction can be expressed as , , , , n X X X X X =  . In this study, F represents the experiment fracture toughness obtained from the X input factors, R denotes the n n × matrix, whose (i, j) entry is  , n is the number of observations sample, and r is the vecto whose ith element is ( ) The Or-K model can be predicting using the unknow functionŷ(X) aŝ where u(X) denotes a global model and ε(X) denotes a local model. The sample of input X is interpolated using a Gaussian random function. The correlation between Z(X i ) and Z X j is related to the distance between the two corresponding points X i and X j . The distance function between point X i and X j is expressed as where θ k 0 ≤ θ k ≤ ∞ is the kth element of the correlation vector parameter θ. The correlation between X i and X j is defined as The Or-K model prediction can be expressed aŝ where . . , f (X n )] T is the output of evaluation function at X = [X 1 , X 2 , X 3 , . . . , X n ]. In this study, F represents the experiment fracture toughness obtained from the X input factors, R denotes the n × n matrix, whose (i, j) entry is Corr Z(X i ), Z X j , n is the number of observations sample, and r is the vector whose ith element is r i (X) = Corr Z(X), Z X i (32) u(X) is assumed to be constant in the Or-K model, andû is given bŷ It is defined as where I denotes an n-dimensional unit vector. The unknown parameter θ in Equation (29) is known as a hyperparameter that affects the prediction performance of the OR-K model, which can be estimated using the maximum likelihood estimation following In this study, the generic algorithm optimization (GA) was used for the optimized maximum likelihood estimation. For given unknows of θ, the variance of Gaussian distribution (σ 2 ) can be defined as where the variables are the same according to the Equation (31).

Brief of Multi-Fidelity Surrogate Model
As mentioned above, for the multi-fidelity surrogate model, data substitution is performed on the basis of the Kriging model. The multi-fidelity surrogate modeling process is shown in Figure 8. The surrogate model for multi-fidelity approaches a radius bias function (RBF) to represent the global model in Equation (37) (Note. High-fidelity data denotes fracture toughness data obtained from experiments and low-fidelity data denotes fracture toughness obtained from the fracture criteria) where local deviations r T R −1 (F h −μ − FR) are evaluated based on the high fidelity data set obtained using the Or- . . , f h (X n )] T and is the output data obtained from high-fidelity function at X = [X 1 , X 2 , X 3 , . . . , X n ] and FR = [ f r (X 1 ), f r (X 2 ), f r (X 3 ), . . . , f r (X n )] T . Note that µ(X) is a mean value of the Gaussian process of the high-fidelity data, assumed to be a constant value expressed by Equation (34), and the definition ofμ is given by Equation (33). The term [µ(X) + f r(X)] is a compound of the Or-K model term of high fidelity data µ(X), which has been defined in Equation (28) and the RBF term of low fidelity data f r (X) predicted from the low fidelity data can be expressed as f r (X) = a 0 + a 1 f l (X) where f l (X) is a function predicted by the RBF using low fidelity data, and a 0 and a 1 are correlation terms between the low fidelity and high fidelity data. The function predicted by the RBF using low fidelity data ( f l (X)) can be expressed as where Φ(X) is an RBF, and w i is a weigh function. A multi-quadratic function is applied as an RBF. The weight is determined from the interpolation conditions where a i,j = Φ x i − x j and F l = [ f l (X 1 ), f l (X 2 ), f l (X 3 ), . . . , f l (X n )] T are the values of the low fidelity function at X = [X 1 , X 2 , X 3 , . . . , X n ].
As mentioned above, for the multi-fidelity surrogate model, data substitution is performed on the basis of the Kriging model. The multi-fidelity surrogate modeling process is shown in Figure 8. The surrogate model for multi-fidelity approaches a radius bias function (RBF) to represent the global model in Equation (37) (Note. High-fidelity data denotes fracture toughness data obtained from experiments and low-fidelity data denotes fracture toughness obtained from the fracture criteria) Figure 8. The modeling process of the multi-fidelity surrogate model. The variance of Gaussian distribution (σ 2 ) can be defined as The unknown parameter (θ, a 0 and a 1 ) of multi fidelity-surrogated model can be estimated using the same maximum likelihood estimation as the Or-K model. According to the various type of low fidelity function (fracture criteria) that were used for the multifidelity surrogate model, the three multi-fidelity models were created. The model namely GM-K represents the multi-fidelity model based on the GMTS criteria, the AS-K model represents multi-fidelity model based on the ASED criteria, and the ER-K model represents the multi-fidelity model based on the GMERR criteria. For the multi-fidelity surrogate model, the data obtained from the experiment fracture toughness were used to create the model along with data obtained from the fracture criteria that was calculated by the separation of thee two thicknesses and mixed mode parameters from 0 to 1.0 with an increment 0.1 for each criterion.

Results and Discussions of the Artificial Intelligence Prediction Models
According to the modeling process of the multi-fidelity surrogate model that combines high-fidelity data, a high precision response to the behavior of the output of the problems and low-fidelity data that show a low precision response to the behavior of the output of problems occur together in the modeling process. When closely considering this modeling method, also called the multiple certainties or multiple precision modeling methods, in which the high-fidelity data (actual experimental fracture toughness) represents the data with high certainty, whereas the low-fidelity data (fracture criteria) is simple to obtain but is inaccurate or has many errors, the data represents a low certainty response to the behavior of the output of the problems. The prediction results of the Or-K model in the modeling process are shown in Figure 9, which shows the modeling based on the experimental fracture toughness. To evaluate all artificial intelligence, the prediction performance metrics in different terms including the percentage-base error measurement (MAPE, Equation (23)), the scale based-error measurement (RMSE, Equation (26)), and the goodness-of-fit statistic (R 2 , Equation (27)) between prediction results of the testing data set obtained via holdout cross-validation method and actual mixed mode I/II fracture toughness results were measured. When considering the prediction performance of the artificial intelligence model in mixed-mode I/II fracture toughness problem, the performances were measured separately into mode I and mode II according to the loading characteristic that causes specimens to fracture. The Or-K model had performance metrics including R 2 , MAPE, and RMSE values equal to 0.896, 16.70%, 0.207 and 0.859, 14.79%, and 0.230 at mode I (K I ) and mode II (K I I ) fracture toughness, respectively. The model performance metrics (MAPE) of the Or-K model demonstrate low prediction performance caused by a fundamental aspect of the artificial intelligence modeling due to the small dataset usage in modeling that affected the learning process (6 data points were used to generate the Or-K model). The prediction results of the multi-fidelity surrogate model, which aims to improve the prediction performance of the artificial intelligence mode in case the data used in the modeling process is limited, is shown in Figures 10-12. The multi-fidelity surrogate model is based on the GMTS criteria ( Figure 10), which used 28 data points in the modeling process (6 data from experiment and 22 data from criteria) and had R 2 , MAPE, and RMSE values at mode I and II fracture toughness equal to 0.954, 12.41%, 0.147, and 0.935, 12.54%, and 0.156, respectively. The multi-fidelity surrogate model is based on the ASED criteria ( Figure 11) and had R 2 , MAPE, and RMSE values at mode I and II fracture toughness equal to 0.933, 15.07%, 0.170 and 0.927, 12.25%, and 0.164, respectively. The multi-fidelity surrogate model is based on the GNERR criteria ( Figure 11) and had R 2 , MAPE, and RMSE values at mode I and II fracture toughness equal to 0.945, 12.23%, 0.155 and 0.961, 12.28%, and 0.122, respectively. When the performance metrics of each multi-fidelity surrogate model were compared, it was found that firstly, for the R 2 values, the GM-K model had the highest mode I fracture toughness and the ER-K model had the highest at mode II fracture toughness. Secondly, for the MAPE values, the ER-K model had the lowest at mode I fracture toughness and the AS-K model had the lowest MAPE at mode II fracture toughness. Finally, for the RMSE values, the GM-K model had the lowest fracture toughness at mode I and the Ri-K model had the lowest fracture toughness at mode II.
The prediction performance metrics of the multi-fidelity surrogate model demonstrate that the accuracy of the model had increased when compared to the Or-K model, which only used data from experiments. The R 2 values of the highest accuracy model had improved when compared to the Or-K of around 6.47% and 11.87%, the MAPE values decrease around 26.77% and 17.17%, and the RMSE values decrease around 28.98% and 46.96% at mode I and II fracture toughness, respectively. When compared to the fracture criteria, unfortunately, the mode I fracture toughness from the multi-fidelity model shows a decrease in predictive performance. The MAPE values of the GM-K, AS-K, and ER-K models were increased by around 15.01%, 46.88%, and 56.39%, respectively, when compared to the based fracture criteria used for each modeling. The error of the multi-fidelity surrogate model at mode I fracture toughness occurs mainly in mixed-mode loading, consistent with the predicted results of the fracture criteria in mixed-mode loading ( Figure 5) used as additional data in modeling, whereas under pure mode I and II loading, fewer errors occur compared with mixed-mode loading because prediction results of the fracture criteria are close to the experimental results because the fracture criteria requires fracture toughness under pure mode loading in the fracture criteria equation. It is for this reason that the multi-fidelity surrogate model has more errors than the fracture criteria; however, when considering the interpretation of the MAPE values, the performance of the three models was also considered to be a good predictor. In mode II fracture toughness, when compared to the fracture criteria, the multi-fidelity surrogated model efficiency increased. When considering the MAPE value, it was found that the MAPE values of the GM-K, AS-K, and ER-K models were decreased around 55.86%, 50.80%, and 48.94%, respectively. The decrease in MAPE showed a marked improvement in predictive performance, when before it was only possible to predict trends of the data when using the fracture criteria.      The prediction performance metrics of the multi-fidelity surrogate model demonstrate that the accuracy of the model had increased when compared to the Or-K model, which only used data from experiments. The 2 R values of the highest accuracy model had improved when compared to the Or-K of around 6.47% and 11.87%, the MAPE values decrease around 26.77% and 17.17%, and the RMSE values decrease around 28.98% and 46.96% at mode I and II fracture toughness, respectively. When compared to the fracture criteria, unfortunately, the mode I fracture toughness from the multi-fidelity model shows a decrease in predictive performance. The MAPE values of the GM-K, AS-K, and ER-K models were increased by around 15.01%, 46.88%, and 56.39%, respectively, when compared to the based fracture criteria used for each modeling. The error of the multi-fidelity surrogate model at mode I fracture toughness occurs mainly in mixedmode loading, consistent with the predicted results of the fracture criteria in mixed-mode loading ( Figure 5) used as additional data in modeling, whereas under pure mode I and II loading, fewer errors occur compared with mixed-mode loading because prediction results of the fracture criteria are close to the experimental results because the fracture criteria requires fracture toughness under pure mode loading in the fracture criteria equation. It is for this reason that the multi-fidelity surrogate model has more errors than the fracture criteria; however, when considering the interpretation of the MAPE values, the performance of the three models was also considered to be a good predictor. In mode II fracture toughness, when compared to the fracture criteria, the multi-fidelity surrogated model efficiency increased. When considering the MAPE value, it was found that the MAPE values of the GM-K, AS-K, and ER-K models were decreased around 55.86%, 50.80%, and 48.94%, respectively. The decrease in MAPE showed a marked improvement in predictive performance, when before it was only possible to predict trends of the data when using the fracture criteria.
The additional fracture toughness was tested to measure the model prediction performance in datasets that are different from the modeling process. The prediction results of the multi-fidelity surrogated model with the same data and different levels of thickness and mode mix parameters are shown in Tables 3 and 4 according to the mode loading. Further experimental results show that the fracture toughness obtained by the prediction model is very close to the experiment results. When considering mode I fracture toughness, the MAPE values of the Or-K model, GM-K model, AS-K model, and ER-K model are equal to 18.90%, 9.42%, 11.35%, and 8.90%, respectively, whereas at mode II, fracture toughness is equal to 9.83%, 3.33%, 3.33%, and 2.80%, respectively. The prediction results in Tables 3 and 4 showed a very close prediction between the Or-K The additional fracture toughness was tested to measure the model prediction performance in datasets that are different from the modeling process. The prediction results of the multi-fidelity surrogated model with the same data and different levels of thickness and mode mix parameters are shown in Tables 3 and 4 according to the mode loading. Further experimental results show that the fracture toughness obtained by the prediction model is very close to the experiment results. When considering mode I fracture toughness, the MAPE values of the Or-K model, GM-K model, AS-K model, and ER-K model are equal to 18.90%, 9.42%, 11.35%, and 8.90%, respectively, whereas at mode II, fracture toughness is equal to 9.83%, 3.33%, 3.33%, and 2.80%, respectively. The prediction results in Tables 3 and 4 showed a very close prediction between the Or-K model and all multifidelity models in the levels of the input factors used to create the model. The prediction results model is very close because of the fact that the model had seen data at this level before in the modeling process. On the other hand, the prediction result of the Or-K model and all multi-fidelity models are quite different in the levels of the input factors that differ in the modeling process, where the multi-fidelity models have values that are clearly close to the experimental results. The above results show that modeling to predict mixed-mode I/II fracture toughness results from combining the experimental data with the fracture criteria to improve model accuracy in case of factor levels never seen before. This indicates that the multi-fidelity surrogated model can solve the prediction issue in the case of predicting fracture toughness of expensive or low-volume materials, making it impossible to perform a large number of tests.

Conclusions
This research aims to improve the prediction performance of artificial intelligence models (Kriging model) on mixed-mode I/II fracture toughness of PMMA, which is modeled using small datasets to respond in cases where data acquisition is limited. Prediction model improvement involves combining data obtained from experiment testing with data obtained from fracture criteria using a multi-fidelity surrogate model. The multi-fidelity surrogate model based on the fracture criteria includes the generalized maximum tangential stress criteria (GMTS), the average strain energy density criteria (ASED), and the generalized maximum energy released rate criteria (GMERR). The results of the research can be summarized as follows:

1.
As for the fracture criteria, the elementary fracture toughness prediction results are very close to the experimental results under pure mode loading since the fracture criteria equations rely on critical fracture toughness under pure load (K IC , K I IC ), which was obtained from the experiment. For the predicted results at mixed-mode loading, the values were found to be rather inconsistent with the experimental results.

2.
As for the original Kriging model, the predicted fracture toughness was rather inaccurate compared to the experimental results in the modeling process. The model had R 2 values equal to 0.896 and 0.859, MAPE values equal to 16.70% and 14.79%, and the RMSE equal to 0.207 and 0.230 when considering modes I and II fracture toughness, respectively. 3.
The prediction performance of the multi-fidelity surrogate model, which is modeled on experimental data as well as the elementary prediction data obtained from the fracture criteria, was found to be higher than that of the original Kriging model or the fracture criteria. For the multi-fidelity surrogate model, the model's performance depends on the fracture criteria used in the modeling process.

4.
The The multi-fidelity surrogate model based on the fracture criteria will ostensibly perform better than the original Kriging model, which solely relied on experimental data in the modeling process, in case the input factors to be predicted differ from the input factors used in the modeling process. 6.
The multi-fidelity models' prediction performance indicated that they are very useful in situations where testing materials are difficult to obtain or prepare for in order to gather enough data to apply artificial intelligence techniques to the problem of fracture toughness. They also assist in reducing the costs associated with data acquisition.