Quantitative Prediction of Fracture Toughness ($K_{{\rm I}c}$) of Polymer by Fractography Using Deep Neural Networks

Fracture surfaces provide various types of information on fracture. For example, a dimple-like pattern indicates the occurrence of plastic deformation during fracture and a consequent requirement of a large amount of energy to fracture the material[1]. However, smooth fracture surfaces as seen in cleavage fractures are suggestive of brittle fracture[1]. From a fracture surface, quantitative parameters of the fracture process, such as the fracture toughness $K_{{\rm I}c}$ or Charpy impact absorption, can be estimated. In previous studies, such parameters have been estimated using the three-dimensional (3D) information of a fracture surface, i.e., its roughness. For example, Mandelbrot et al. obtained the 3D landscape of a fracture surface by observing its structural changes while continuously slicing the sample. They found a correspondence between the Charpy impact absorption and the fractal dimension by the fractal analysis of this contour structure[2]. Another study revealed that $K_{{\rm I}c}$ can be estimated by extracting the roughness of a fracture surface measured by stereo imaging using an electron microscope[2, 3]. From an engineering perspective, it is time-consuming and expensive to obtain the 3D information of a fracture surface; thus, it is desirable to estimate parameters of fracture process such as $K_{{\rm I}c}$ from two-dimensional (2D) images, which can be acquired at a relatively low cost. The results of some previous studies have suggested that the contrast in 2D images contains 3D roughness information and that the identification of the fracture type, such as ductile or brittle fracture, is possible[1, 4, 5]. However, this contrast is also affected by the angle of the measurement target, the type of observation probe (such as light or electrons), and the interaction between the probe and the measured sample. Therefore, in general, extracting a 3D structure from contrast is difficult in principle. Thus, compared with previous methods that rely on 3D information, estimating parameters of fracture process such as $K_{{\rm I}c}$ from 2D images is difficult.

In this study, we establish a framework for fractography by which we can directly estimate $K_{{\rm I}c}$ from 2D fracture surface images using a deep neural network (DNN) by the transfer learning (TL) method. In recent years, deep learning methods to estimate the depth position of an object from a single 2D image have been rapidly developed[6, 7]. Thus, deep learning can extract the features of an uneven structure from the contrast of its 2D images under the effects of the above-mentioned complex factors. These recent advances in deep learning research suggest that $K_{{\rm I}c}$ can be directly estimated from a 2D fracture surface image by extracting its features through deep learning. Typically, image recognition using a DNN requires a tremendous amount of image data, but it is difficult to acquire such a large amount of image data for fractography owing to the high cost of experimentation involving specimen preparation, observation, and the measurement of fracture properties. The dataset utilized in this study contains 770 2D fracture surface images with corresponding $K_{{\rm I}c}$; this is a large number for such types of data but quite small as a training dataset for a regular DNN, as subsequently described in detail. To compensate for the limited data, we used the TL method, which allowed us to construct high-performance prediction models even with a small dataset by transferring machine learning models trained using other large datasets. TL can be used to identify fracture surfaces using DNNs[8, 9, 10, 11, 12]. Thus, it is expected that by using a DNN and TL, a framework for directly estimating $K_{{\rm I}c}$ can be established. Concretely, in this study, we examined whether a DNN can extract the features required for predicting $K_{{\rm I}c}$ from images of fracture surfaces of polymers. The present fracture dataset of polymers, which was collected by the National Institute of Technology and Evaluation[13], includes 770 data items, each being a set of one fracture surface image and its $K_{{\rm I}c}$. This dataset contains macroscale fracture surface data generated by various fracture processes, such as brittle and ductile fracture processes. The results obtained in this study showed that even with such limited data, $K_{{\rm I}c}$ can be estimated from 2D fracture surface images using the established DNN and TL framework.

The remainder of this paper is organized as follows. First, in Sec. II, we describe the fracture toughness test data of the polymeric materials used in this study. In Sec. III, the DNN model used for TL and the regression model are explained. In Sec. IV, the results of the analysis are presented, and in Sec. V, the results are discussed and summarized.

In this study, the “Polymer Fracture Database”[13] is used, which is provided by the National Institute of Technology and Evaluation (NITE). The development of this database as a collaborative effort of NITE, Yamagata University, Osaka Municipal Technical Research Institute, and Meiji University as a project commissioned by the New Energy and Industrial Technology Development Organization (NEDO) began in 1996. Since 2020, the database has been available on the Materials Data Repositry (MDR) managed by the National Institute for Materials Science (NIMS).

This database was established to develop a widely used technological infrastructure to enhance the safety of materials, such as the development of technologies for investigating the causes of accidents arising from the fracture or failure of plastic products. The database contains test data for various polymeric materials used in society, which were obtained from the fracture toughness, tensile, fatigue, creep, essential work of fracture (EWF), and scratch tests. In this paper, data from the fracture toughness test were used, because such data (upper figure of Fig. 1) are the most comprehensive. For the fracture toughness test, data were obtained for various types of samples with different thicknesses, compositions, grades, molding methods, mold temperatures, injection speeds, cylinder temperatures, and with/without welds, at different test temperatures and speeds. The test temperatures were centered at room temperature and generally ranged from $-40$ °C to $+60$ °C. The database contains macro-photographs taken of about a 10 mm square of the fracture surface for each set of the above conditions (lower figures of Fig. 1). One fracture surface image was included for each of the 1248 test conditions. Among the parameters in the database that are expected to be estimable, as mentioned in Introduction, $K_{{\rm I}c}$ for plane strain was selected in this study because it is the main parameter that one would aim to obtain in a fracture toughness test, and it has always been recorded in many tests in the database.

The database includes comma separated values (CSV) files listing information for all tests, PDF files containing fracture surface photos for each set of test conditions, and tab separated values (TSV) files containing load-displacement curve data. Each test was assigned an identification code. Since the fracture surface images were recorded in the lossy compressed JPEG format in the PDF file of each test and bitmap data without degradation could not be obtained, the image data were extracted, as it is in the JPEG format, using the pdfimages[14] command. Note that the extracted image file had its left and right sides reversed from when the PDF was displayed, so the left and right sides were reversed using the convert command of Imagemagick[15]. From among all test data, 770 data records were obtained by extracting test data in which $K_{{\rm I}c}$ was recorded and which included fracture surface data. These 770 data records were used as the dataset for the analysis in this study.

In this section, we describe the method adopted for building a predictive model of $K_{{\rm I}c}$ using feature-extraction-type transfer learning[16], which is a combination of a DNN trained on natural images and Gaussian process regression[17]. Specifically, we developed a model that uses a DNN as a feature extractor and regresses its features and $K_{{\rm I}c}$ by Gaussian process regression.

The DNN model used in this study was the vgg16 model[18], which identifies natural images such as animals and vehicles. vgg16 is composed of convolution, pooling, and fully connected layers. A convolution layer convolves the input space using a convolution filter. Convolution filters of DNNs trained on natural images frequently behave similarly to Gabor filters. A Gabor filter can extract image features such as the edges of graphic structures in images. In vgg16, a structure comprising a stack of convolution layers alternated with pooling layers is formed, and the pooling layers compress the dimensions by replacing certain regions of the convolution layers with representative values such as the maximum and average of the region. Pooling can be viewed as a transformation that renders the input space coarse grained. Thus, after training on natural images, vgg16 is expected to extract useful features on various scales for natural images. In the prediction of $K_{{\rm I}c}$ using fracture surface images, the evaluation of the edges of the graphic structures in the images is also considered important. Depth estimation research using images has shown that using a DNN pretrained with the natural image dataset called ImageNet[19] as a feature extractor improves the accuracy of depth estimation after fine-tuning, compared with the method without pretraining[7]. This finding suggests that DNNs trained on Imagenet can extract the image features necessary for depth estimation. As noted in Introduction, information on depth is useful in predicting $K_{{\rm I}c}$. Thus, the image features extracted by vgg16 trained on natural images at each convolution layer are considered to yield useful features for predicting $K_{{\rm I}c}$. Therefore, we employed feature-extraction-type transfer learning[16], that used vgg16 trained on natural images as an extractor of image features to construct a prediction model of $K_{{\rm I}c}$. With such transfer learning, a framework can be constructed for learning a large amount of data at a lower cost owing to the exclusion of the fine-tuning of the DNN. However, in the absence of fine-tuning, if certain image features required for the prediction of $K_{{\rm I}c}$ are not included in the set of image features required for the classification of Imagenet, then they remain missing. In this study, we examined whether the prediction of $K_{{\rm I}c}$ is still possible under such conditions.

The predictive model of $K_{{\rm I}c}$ was constructed as follows. First, we prepared vgg16 (top figure of Fig. 2) trained on the ImageNet dataset. In this study, publicly available learned weights[20] were used as the trained vgg16. In addition, a library[21] using TensorFlow was adopted as a prototype for the vgg16 implementation. When an image is input to the pretrained vgg16, the features are obtained as the outputs of the five pooling and three fully connected layers referred to as Layer 1–Layer 5 and Fc 1–Fc 3, respectively. We denote one of these DNN features as $m\in\{{\rm Layer\>1,Layer\>2,Layer\>3,Layer\>4,Layer\>5,Fc\>1,Fc\>2,Fc\>3}\}$. Let $\mathbf{h}^{i}_{m}$ be a $d_{m}$-dimensional DNN feature for a certain fracture image $i$. Because the number of data $N=770$ was much smaller than $d_{m}$, $\mathbf{h}^{i}_{m}$ was compressed to $N$ dimensions by principal component analysis (PCA) to reduce the computational complexity. Let $\mathbf{x}^{i}_{m}$ be a compressed DNN feature. Note that in a PCA that compresses data into the dimension of the sample size, information on the covariance of the dataset is not lost. Therefore, a regression model should be learned to predict $K_{{\rm I}c}$ perfectly consistent with that in the uncompressed case in Gaussian process regression. Each fracture surface image, which was approximately $400\times 400$ pixels in size, was prescaled using zero-order spline interpolation[22] to ensure that it had the same number of pixels as the natural image ($224\times 224$). A predictive model of $K_{{\rm I}c}$ was constructed by Gaussian process regression using the dimension-compressed DNN feature $\mathbf{x}^{i}_{m}$ as the explanatory variable (bottom figure of Fig. 2). Note that in the entire regression model consisting of vgg16 and Gaussian process regression, only the regression parameters of the latter were trained on the fracture surface image data, and the DNN model used the parameters from the pretraining as is. Gaussian process regression allows the flexible construction of regression models by appropriately choosing a kernel $k(\cdot,\cdot)$. In this study, the rational quadratic kernel

and the Gaussian kernel

were the candidates for $k(\cdot,\cdot)$. Here, $\theta_{k}$, which determines the shape of the kernel, is the hyperparameters of the Gaussian process regression.

To construct such a regression model, the DNN feature $m$, the kernel function type $k$, and the kernel function hyperparameters $\theta_{k}$ must be properly selected. In this study, these parameters are selected on the basis of a Bayesian inference framework[23]. For the selection of $m$ and $k$, the Bayesian model selection framework is employed. On the basis of an indicator called Bayesian free energy[24], the Bayesian model selection framework selects the model that considers discrete-valued hyperparameters such as $m$ or $k$ (see Appendix A). In general, when using machine learning to build predictive models, it is necessary to prevent over-fitting, where the model over-fits the training data and loses predictive performance. Using Bayesian free energy as an indicator for model selection, it is possible to select a model in which over-fitting does not occur. $\theta_{k}$, that is, the continuous-valued hyperparameter, can also be estimated on the basis of Bayesian free energy (see Appendix A). Thus, the model and hyperparameters are selected so as to have a lower Bayesian free energy. $\theta_{k}$ was optimized using the L-BFGS-B algorithm [25], a type of quasi-Newtonian method. Newton’s method sometimes gives a local solution depending on the initial values of the parameters to be optimized. Therefore, in this study, we searched for initial values that minimize the Bayesian free energy in the range shown in Table 1. The brackets $(\>,\>)$ to the right of each parameter in Table 1 indicate how to set the grid. The first element in the brackets is a scale of the space to be gridded and the second element is the number of grids.

For comparison with the case where the DNN feature is not used, a regression of $K_{{\rm I}c}$ was performed directly from the fracture surface images by Gaussian process regression with the same procedure as in the case of DNN features. Hereafter, the obtained regression model is referred to as the “direct regression” model.

As a result of Bayesian model selection, Layer 3 with a rational quadratic kernel is the optimal model for predicting $K_{{\rm I}c}$, where the Bayesian free energy is minimum (Fig. 3). For all DNN features, the Bayesian free energy with the rational quadratic kernel is always lower than that with the Gaussian kernel. This indicates that the rational quadratic kernel is always a plausible kernel type for any DNN feature. In terms of optimal DNN features, Layers 3 and 4 were confirmed to be plausible features specifically with a lower Bayesian free energy than the others, regardless of the kernel type. The most plausible regression model with the minimum Bayesian free energy was obtained when the rational quadratic kernel was selected as the kernel function and Layer 3 for the DNN features. The difference in the Bayesian free energy of the regression model using Layer 3 and Layer 4 features in the rational quadratic kernel was only about 3.36. However, since Bayesian free energy is defined as the logarithm of the likelihood function, which expresses the likelihood of the model (Appendix A), the features of Layer 3 are more than $\exp(3.36)=28$ times more likely than those of Layer 4 when compared in the original linear space. Therefore, a model that performs regression using a rational quadratic kernel with Layer 3 is selected here. This selected regression model has a lower Bayesian free energy than direct regression. It is noteworthy that under the same kernel function, there are DNN features that yield regression models with lower performance in terms of Bayesian free energy than direct regression. For example, in the rational quadratic kernel, only Layers 3 and 4 had a lower Bayesian free energy than direct regression, and other DNN features such as convolution and fc layers produced regression models with a higher Bayesian free energy than direct regression. Under a feature extraction type of transfer learning, higher-layer features than the last layer of the convolution layer[16, 26, 27, 28, 29], Layer 5 in our model, are often used, but in the case of this study, they were not considered. In other words, depending on the material to be analyzed and the mechanical properties to be estimated, transfer learning would not be effective unless the appropriate Layer is selected. The selected Layer 3 corresponds to a resolution of about 10$\%$ of the input image. This means that the spatial scale of the best DNN feature is on a relatively coarse-grained scale. It is suggested that a complex feature, different from the simple Gabor filter, is useful for $K_{{\rm I}c}$ estimation.

The regression of $K_{{\rm I}c}$ was performed using the DNN features, the kernel type, and the kernel hyperparameters $\theta_{k}$ selected in terms of Bayesian free energy. Figure 4 shows the regression results. To use all the data effectively, we estimated the prediction performance on the basis of the framework of the 20-fold cross-validation method. The black dots are the integrated prediction results for test data in all test folds where data are not used for training. The horizontal axis is the measured value and the vertical axis is the predicted value of the constructed model. The left figure shows the result of predicting $K_{{\rm I}c}$ directly from the fracture image (direct regression model). The direct regression returns a constant of around 3 for most of the cases, and there is virtually no correlation between the predicted and measured values. This indicates that the direct regression method cannot extract enough information to predict $K_{{\rm I}c}$. Note that the regression model using DNN feature extraction with Layer 3 has high accuracy, as shown by Fig. 4(b). The $R^{2}$ value of the prediction model using Layer 3 features was 0.65 (correlation coefficient is 0.81). We also found that our model can predict $K_{{\rm I}c}$ with a standard deviation of the estimation error of about $\pm$0.37 [MPa$\sqrt{m}$] in the range of around 1–5 [MPa$\sqrt{m}$] (Fig. 4). The high accuracy demonstrates that it is possible to predict $K_{{\rm I}c}$ using only the 2D fracture image without 3D height information. The high accuracy also suggests that the morphology of the fracture surface contains some information about resistance to fracture.

Thus, the comparison with the poor result of the direct prediction method proves that the transfer learning using a DNN trained on natural images as a feature extractor is crucial in building a predictor for $K_{{\rm I}c}$ that has high estimation accuracy. It was also revealed that the importance of the proper use of DNN information, i.e., the selective use of Layer 3 in this case, should be emphasized.

The present result demonstrates that the DNN trained with TL opens the way to quantitative fractography in which one can estimate the parameters of fracture process from the fracture surface image even with a small dataset. In this study, the DNN pretrained on the Imagenet dataset was used directly as a feature extractor. Fine-tuning the DNN with fracture surface image data is expected to improve the performance of the model obtained in this study. It is also assumed that a higher estimation accuracy of $K_{{\rm I}c}$ could be achieved by using not only the fracture surface image but also the additional explanatory variables describing the target material, such as the production conditions or the operating environment.

The proposed framework, which enables the construction of regression models in a few hours, can be used for the low-cost screening of profitable candidates on which to apply machine learning analysis to estimate values representing mechanical property from a large number of image datasets available in the field of materials science. Therefore, if the utilization of the proposed framework is expanded, it is expected that machine learning methods in materials science will be used more broadly, quickly, and efficiently.

This work was supported by the Council for Science, Technology and Innovation (CSTI), Cross-Ministerial Strategic Innovation Promotion Program (SIP), “Structural Materials for Innovation” and “Materials Integration for Revolutionary Design System of Structural Materials” (funding agency: JST), JST PRESTO Grant Number JPMJPR212A, and JSPS KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas “Discrete Geometry for Exploring Next-Generation Materials” 20H04648 and 22K13979.

Let $\mathbf{x}^{i}_{m}$ be the N-dimensional compressed output of the $m$-th layer using PCA when a certain fracture surface image $i$ is input to the pretrained DNN. Also, let the $K_{{\rm I}c}$ value of fracture surface image $i$ be $y^{i}$. Suppose there are $N$ samples of $\mathbf{x}^{i}_{m}$ and $y^{i}$. Let them form the matrix $\mathbf{X}:=(\mathbf{x}^{1}_{m},\mathbf{x}^{2}_{m}\cdots\mathbf{x}^{N}_{m})$ and the vector $\mathbf{Y}:=(y_{1},y_{2}\cdots y_{N})^{T}$.

Unlike regression based on the maximum likelihood method, Gaussian process regression does not involve the estimation of parameters of the regression function, such as regression coefficients in linear regression. In Gaussian process regression, the regression function $\mathbf{F}:=\mathbf{F}(\mathbf{X})=(f(\mathbf{x}^{1}_{m}),f(\mathbf{x}^{2}_{m})\cdots f(\mathbf{x}^{N}_{m}))$ itself is the target of estimation, and a likelihood function in which $\mathbf{F}$ is the probability variable is introduced as

where $Z$ is the normalization constant and $\mathbf{I}$ is the identity matrix. From Bayes’ theorem $P(B|A)=\frac{P(A|B)P(B)}{P(A)}$, the regression function $\mathbf{F}$ is estimated as the posterior probability distribution ${\rm P}(\mathbf{F}|\mathbf{Y})$ as

where the prior distribution ${\rm P}(\mathbf{F})$ is assumed to be an $N$-dimensional normal distribution with mean 0 and variance–covariance matrix $\mathbf{K}$, and $k(\mathbf{x}^{i}_{m},\mathbf{x}^{j}_{m})$ is a kernel function that indicates the similarity between $\mathbf{x}^{i}_{m}$ and $\mathbf{x}^{j}_{m}$. In this study, the rational quadratic kernel

and the Gaussian kernel

are employed as candidates for the kernel function. The hyperparameters of each type of kernel are defined as $\theta_{\rm Rational}=\{\theta_{0},\theta_{1},\theta_{2}\}$ and $\theta_{\rm Gauss}=\{\theta_{0},\theta_{1}\}$. By reorganizing Eq. (6), it is revealed that the posterior distribution is a normal distribution with respect to $\mathbf{F}$.

where $Z^{\prime}$ is the normalization constant. This type of regression in which the prior and posterior distributions are normal distributions is referred to as Gaussian process regression. Using the obtained regression model ${\rm P}(\mathbf{F}|\mathbf{Y})$ [Eq. (9)], the predicted value $y^{\rm new}$ at the unknown point $\mathbf{x}_{m}^{\rm new}$ is given as the probability ${\rm P}(y^{\rm new}|\mathbf{Y},m,k,\theta_{k})$:

where the prior distribution ${\rm P}(m,k)$ is assumed to be uniformly distributed.

In this study, on the basis of the framework of Bayesian model selection and empirical Bayesian methods[23], the DNN feature $m$, the kernel function type $k$, and its hyperparameters $\theta_{\rm k}$ were determined by maximizing the following marginal likelihood functions called Bayesian evidence:

In the Bayesian model selection framework, the likelihood of a model or hyper-parameter is often compared not with the marginal likelihood itself, but with its negative logarithm called Bayesian free energy[24],

where $k$ denotes the kernel type. This Bayesian free energy is also employed in this study.

Since $\theta_{\rm k}$ of the kernel function has continuous values, it is not possible to search for the value that minimizes the Bayesian free energy in a grid search like $m$ or $k$. Therefore, in this study, the L-BFGS-B algorithm [25], a kind of quasi-Newtonian method, is used to find $\theta_{\rm k}$ that minimizes the Bayesian free energy. Newton’s method becomes a local solution depending on the initial values of $\theta_{\rm k}$. In this study, we explored the initial values that minimize $F(m,k,\theta_{\rm k})$ in the range shown in Table 1. The model parameters $k$ and $m$ were determined with $\theta_{\rm k}$ optimized in this manner.