# Fileset

[Size-dependency_Ito_et_al_R1_unmarked.docx](https://mdr.nims.go.jp/filesets/7aa775ef-87a2-43e8-beff-a7a51fadcfd9/download)

## Creator

Chihiro Ito, Taiyo Maeda, Ryunosuke Higashi, [Toshio Osada](https://orcid.org/0000-0003-1539-9264), Takuma Kohata, Shingo Ozaki

## Rights

[Creative Commons BY-NC-ND Attribution-NonCommercial-NoDerivs 4.0 International](https://creativecommons.org/licenses/by-nc-nd/4.0/)

## Other metadata

[Application of extreme value statistics to internal pore distribution in ceramics and prediction of size dependency of strength scatter](https://mdr.nims.go.jp/datasets/f98d7067-3c3f-4e54-8ce9-407bc59ec36a)

## Fulltext

Application of extreme value statistics to internal pore distribution in ceramics and prediction of size dependency of strength scatterChihiro Itoa,+, Taiyo Maedaa,+, Ryunosuke Higashia, Toshio Osadab,c, Takuma Kohatab, Shingo Ozakib,c,*aGraduate School of Engineering Science, Yokohama National University, Tokiwadai 79-5, Hodogaya-ku, Yokohama 240-8501, JapanbHigh Temperature Materials Group, Research Center for Structural Materials, National Institute for Materials Science, Sengen 1-2-1 Tsukuba, Ibaraki 305-0047, JapancDivision of System Research, Faculty of Engineering, Yokohama National University, Tokiwadai 79-5, Hodogaya-ku, Yokohama 240-8501, Japanc, * Correspondence: E-mail address: s-ozaki@ynu.ac.jp, Tel/Fax No.: +81-45-339-3881+Co-first authors: These two authors contributed equally to this work.AbstractTo implement ceramics in high-reliability components, it is necessary to understand the strength scatter caused by the stochastic distribution of defects in ceramics. In this study, we propose a numerical simulation method to predict the strength scatter of ceramics and its size dependency using microstructural data obtained using X-ray computed tomography and scanning electron microscopy. The generalized Pareto (GP) model, a type of extreme value statistic, was applied to treat the three-dimensional distribution of the pore size. Further, we predicted the bending strengths of the four types of bending tests with different effective volumes, and the results corresponded with the experimental results. Additionally, the method was applied to specimen models with different discretization sizes to verify the advantages of the GP model. The results show the effectiveness of the proposed method for the specimen size dependence of the bending strength and specimen discretization size independence.Keywords: Strength analysis, defects, extreme value statistics, Weibull distribution, size dependency1. IntroductionCeramics have a higher specific strength than general metallic materials and exhibit several superior properties such as heat, acid, and wear resistance. Recently, they have been used for different applications requiring high reliability, such as electronic components (multilayer capacitors)[1–4], biomaterials (dental implants)[5–7], and structural components (aircraft engines) [8–11]. Generally, the fracture toughness of ceramics is lower than that of metallic materials, and ceramics exhibit a stochastic brittle fracture behavior depending on the distribution of intrinsic defects[12–14]. Hence, the Weibull distribution (or plot) of fracture strength has been widely employed to discuss the strength variation in ceramics[15–22].The Weibull distribution is useful for examining the fracture probability of the target specimens or components. However, the fracture strength of ceramics depends on the size and shape of the material and the loading conditions[18,23]. Although it is possible to estimate the size dependence of the fracture strength based on Weibull statistics and effective volume, the lower limit of the fracture strength or safety factor is difficult to determine because the Weibull modulus (shape parameter) is neither material-specific nor accompanied by upper/lower limits and confidence intervals[24].The complementary use of standardized test results, microstructural observations, and numerical simulations to predict the strengths of specimens of different sizes and shapes under arbitrary working load conditions would help solve the abovementioned problems. Hence, a method has been reported for predicting the strength variation (Weibull distribution) of specimens by setting up the spatial distribution characteristics of defects, which are possible fracture initiation points, and performing Monte Carlo simulations[17,18,20,25–27]. Ozaki et al. proposed a method for predicting the strength scatter of ceramics within the framework of finite element analysis (FEA) by incorporating microstructural features into the parameters of a constitutive equation using a fracture mechanics model[28–30]. This permitted the variation in bulk fracture strength to be discussed in relation to microstructural features such as pores and grain boundaries. Takeo et al. applied the FEA method to five rectangular specimens of different sizes and suggested its possible application in predicting the size dependence of the fracture strength of ceramics[31].Furthermore, Ozaki et al. proposed an extended method to utilize information on the defect density (number of pores per unit volume) obtained using scanning electron microscopy (SEM) and X-ray computed tomography (CT)[32]. However, there exist numerous types of pores of various sizes in ceramics, as shown in Figure 1(a). Hence, the method for handling the probability of pore existence must be refined. Additionally, different types of pores exist depending on the generation mechanism, such as pores within grains and on boundaries, medium-sized pores on triple points, and large defects comprising agglomerated pores around the unsintered powders, as illustrated in Figure 1(b)[33,34]. Hence, treating pores of different origins as belonging to the same population is not appropriate, thereby requiring a more sophisticated method. In this study, we used extreme value statistics as a statistical method to focus only on pores on the extreme side of the probability density function (Figure 1(a)), which could be fracture initiation points. The application of extreme value statistics has been widely studied for metallic materials. Numerous attempts have been made to estimate scatter in fatigue strengths[35–39]. Beretta[38] and Beretta et al.[39] summarized the basic concepts for treating the distribution of inclusions in metals using the generalized extreme value (GEV) model for the block maxima concept and generalized Pareto (GP) model for the peak-over threshold concept, respectively. Further, methods for dealing with the presence of different types of inclusions are discussed. Meanwhile, the application of extreme value statistics to ceramic strength scatter has also been attempted[27,40,41]. Chao et al. used the GEV distribution of the estimated pore size from the measured pore size distribution. They reported that the estimated pore size could be combined with linear fracture mechanics to predict the fracture strength distribution of sintered silicon nitride[40].This study aims to develop a numerical simulation method for the strength scattering of ceramics of different sizes and shapes based on the results of microstructural observations. Hence, we propose a universal method that resolves the dependence not only on the size and shape of the component but also on its discretization. The volume of the discretized region (element) was considered as one block, which was used as the basis for the extreme value statistics. The GP model[42], a type of extreme value statistic, was applied to describe the number density and probability density function for the size of only those pores that may be candidates for the fracture initiation point.The remainder of this paper is organized as follows: Sections 2 and 3 present the strength and microstructural analyses, respectively, Section 4 describes the proposed simulation method, and Section 5 presents the simulation results and discusses the validity of the proposed method, including its predictive performance for the lower strength limit. Section 6 presents the conclusions of this study.2. Materials and strength testIn this study, high-purity alumina AS999 (Ferrotec Material Technologies Corporation, Japan) was used as the target material for experiments and simulations. To determine the size dependence of fracture strength, four types of bending tests with different effective volumes were performed by varying the external and internal span lengths. These comprise two types of three-point bending tests (3pb-S: external span length of 16 mm and 3pb-L: external span length of 30 mm) and two types of four-point bending tests (4pb-S: external span length of 30 mm and 4pb-L: external span length of 60 mm) with internal span lengths of 10 mm. The test method followed the Japanese Industrial Standard (JIS R 1601), and the surfaces of the specimens were mirror-polished. The geometries of the specimens are presented in Table 1. All the specimens used for the bending tests were obtained from the same lot.Bending tests were performed at a crosshead speed of 0.5 mm/s at room temperature. The bending strength  was calculated using the following equation:   where P is the peak value of the jig reaction force, S1 the external span length, S2 the internal span length, t the specimen width, and W the specimen thickness, as shown in Figure 2(a).The bending strength of each specimen was determined using the two-parameter Weibull distribution:  where m is the Weibull modulus (m-value) and  the scale parameter (-value). We adopted the median rank method to calculate the cumulative probability . The number of specimens was  for each size condition.The following equation was used to organize the results based on the effective volume [43]:  where  is the volume of the specimen wherein bending stress is generated. In this study, we evaluated  by substituting the experimental or calculated m-values for each bending test geometry with S1, S2, t, and W into Eq. .For the size dependence of the -value, the following equation is used if the m-value is independent of [17]:  where  denotes the scale parameter when the effective volume is V0. In this study, we evaluated the size dependence using Eq.  for comparison with the experimental and simulation results. Although the experimental m-value varied depending on , the m-value obtained from 3pb-S, which had the smallest  of all the tests, was used in Eq. .The fracture toughness KIC of the target material is required to predict the bending strength based on the fracture mechanics model, as shown in Figure 2(b). The fracture toughness was evaluated using a previously reported method[44,45]and determined by the nonlinear relationship between the three-point bending strength of the cracked specimens and size of a semicircular surface crack ranging from 20 to 310 μm. This was introduced by Vickers indentation with different indent forces ranging from 9.81 to 292.2 N. The fracture toughness of AS999 used in this study was determined to be . 3. Microstructural feature analysisThis section describes the measurement method for the microstructural data (grain size, pore size, and pore aspect ratio) and the approximation using the respective probability density functions, which serve as input data for the numerical simulation of the fracture strength.3.1 2D/3D microstructure observationsFor detailed pore type identification, the microstructure within the wide region on XE −YE plane parallel to both the specimen longitudinal and thickness directions shown in Figure 2(a) was observed using a scanning electron microscope (SEM; ZEISS Gemini SEM 300, Germany) equipped with an electron backscattered diffraction (EBSD) pattern method. As shown in Figures 3(a) and (b), the cross-sectional area (4.085×3.08 mm) was divided into 18×18 sections, and 324 high-resolution SEM images were automatically captured on a bending specimen. The image processing software ImageJ was used to process the pore features, which were classified into two types: “small pores within grains and on boundaries” (hereinafter referred to as small pores) and “medium pores on triple points + large defects comprising agglomerated pores around unsintered powders” (hereinafter referred to as large pores), as shown in Figure 3(c). The machine-learning function (Trainable Weka Segmentation method) in ImageJ software was used for classification. Data with small-sized pores within grains and on grain boundaries can be detected owing to the high resolution of FE-SEM. Furthermore, the two-dimensional (2D) distribution characteristics of the major pore radii R and pore aspect ratio A (ratio of the minor radius to the major radius) were analyzed using an elliptical approximation. The grain size distribution was obtained using the TSL OIM data collection program after the EBSD analysis. The distribution of the local relative density of AS999 was narrow and it was treated as a constant value in this study.Figure 4 shows the fracture initiation sites in the lowest-strength specimen (bending strength of 335 MPa) in the three-point bending tests (3pb-S). The fracture initiation sites of all the test specimens were confirmed to be internal defects, which corresponds to the large pores (Figure 3(c)).To obtain the three-dimensional (3D) spatial distribution of the pores in the specimen, direct observations were performed using high-resolution X-ray CT (ZEISS Xradia 520 Versa), as shown in Figure 5. The observation volume was 279.85×214.02×220.52 µm. The obtained data were evaluated using the analysis software VGStudio Max 3.3 to obtain the 3D distribution characteristics. A total of 1399 pores (min: 3.99 µm, max: 44.33 µm) were measured using X-Ray CT.Owing to the limitations of the measurement volume and resolution of the X-ray CT, it is difficult to assess whether the 3D distributional characteristics of the population are adequately represented using only CT data. Thus, SEM observations with high-resolution and large-area are valuable to obtain distributional characteristics[33,34]. Hence, the 3D distribution data of the pores were evaluated using the SEM images to complement the results of the X-ray CT measurements.Figure 6(a) shows the 2D size distributions of each type of pore (small and large pores) and their combination, which is classified using the ImageJ software (see Figure 3(c)). As schematically illustrated in Figure 1(a), the distribution characteristics differ depending on the generation mechanism. The size distribution of all the pores created without distinguishing between the types had a peculiar shape. Hence, it is necessary to use data only on the pores, which may be a possible fracture initiation point, to properly approximate the probability density function of the size and determine the number density. We converted the 2D distribution to 3D using the Schwartz-Saltykov (SS) method, a type of stereology[46,47] wherein the number of pores in the size group j per unit volume NV is expressed by the following equation using the number of pores per unit area NA:  where k is the total number of size groups and .  is the maximum major pore radius in the observed area, which was 44.33 µm in this study. The following equation for coefficient  was proposed by Takahashi and Suito[46]:  Generally, the accuracy of the converted 3D distribution improves with the total number of size groups k when sufficient data are available[46]. Hence,  was adopted.Figure 6(b) shows the 3D distribution of the pore size (major pore radius R) per unit volume converted using the SS method. The 3D distribution of the pore size per unit volume obtained by X-ray CT is also shown in the figure. The pore data obtained using X-ray CT represents the distribution characteristics of relatively large major pore radii in ceramics which shows significantly fewer pores below 5 µm. This is owing to the removal of pores below 50 µm3 (3.68 m) corresponding to the lower instrumental resolution to reduce the processing costs. The obtained size distribution showed a decrease in the pore number with an increase in the pore size. It was confirmed that the distribution of the SEM data corresponded with that obtained using X-ray CT in the large pore region.3.2 Approximation using probability density functions3.2.1　Grain size Figure 7(a) shows an inverse pole figure of AS999 obtained via EBSD. Similar to previous studies[28–32,48], the distribution of the grain size c arranged by the area ratio was approximated by a log-normal distribution. The best-fit results obtained using cumulative distribution are shown in Figure 7(b). 3.2.2　Pore sizeThe distribution of the major pore radius R was approximated using the GP model. Here, we examine the applicability of the GP model when pores of different generation mechanisms are mixed based on the data obtained using X-ray CT because 3D number density information for large pores is more accurate.The model estimated extreme values using the GP distribution for the right hem (large side) of the population. The cumulative distribution function H and probability density function h of the GP distribution are expressed using the following equations[49]:    In the aforementioned equations, u represents the threshold, and  and  the parameters of the GP distribution, which represent the scale and shape of the distribution function, respectively.To fit the 3D distribution of the major pore radii to the GP distribution, the threshold u was first determined using the sample mean excess plot, which is expressed as follows:  where nu is the number of data points above the threshold, Ri the major pore radius, and Rmax the maximum value of the major pore radius. When fitting the GP distribution, the sample mean excess plot of the random variable Ri varies linearly above the threshold u. The larger the amount of data used, the higher the accuracy. In this study, the threshold was set to the smallest u for which the sample mean excess plot could be approximated using a straight line.Figure 8 shows a sample mean excess plot of the major pore radii obtained from the X-ray CT results. Here, the slope of the plot changes as the threshold increases above 5.9 µm, while the plot can be approximated linearly. The slope further changes above 14 µm. However, this is owing to the limited pores having this measured size, which is insufficient to create the sample mean excess plot. Hence, we set the threshold u to 5.9 µm. The number density of the pores above the threshold was obtained using . This threshold excludes pores that are irrelevant to the fracture initiation point, as shown in Figure 6. This indicates that the GP model employed in this study can extract only the data on candidates of fracture initiation points categorized as “medium pores on triple point” and “large defect comprising agglomerated pores around unsintered powders.” A detailed fit of the GP model to the pore-size distribution of AS999 is provided in the Appendix.The probability weight moment (PWM) method was used to fit the X-ray CT measurement results to a GP distribution[49]. Figures 9(a) and (b) show the fitting results of the major pore radii in the probability density and exponential probability papers, respectively. The GP distribution estimated using the PWM method reproduced the 3D distribution of the pores obtained using X-ray CT (pore distribution A). However, it is important to consider the presence of large pores that cause strength reduction when numerically predicting the scatter of the ceramic strength and its lower limit. Hence, another distribution was prepared in addition to distribution A by adjusting the scaling parameter σ to match the upper side of the distribution in the probability density (pore distribution B). Furthermore, a distribution caught in between distributions A and B was prepared for verification (pore distribution C).3.2.3　Pore aspect ratioThe pore aspect ratio A was estimated using the 2D SEM data. The shapes of the pores were approximated as elliptical for this process. Referring to the threshold obtained when approximating the GP distribution of the major radius, only the pores with a major radius R greater than 5.9 µm were extracted and fitted. Similar to previous studies[28–32], a normal distribution was adopted. The best-fit results obtained using the cumulative distribution are shown in Figure 10.4. Analysis model of bending testsIn this study, the general-purpose interpreter programming language, Python, was used to simulate the bending tests. First, the specimens were discretized using appropriate cubic element sizes to evaluate the bending stress distribution at each location during the bending test, as shown in Figure 2(a). The bending stress of each element, , inside the inner span during a four-point bending test is expressed as  where PR denotes the jig reaction force, I the moment of inertia of the area, and YE the distance from the neutral plane. The bending stress in an element outside the inner span or during a three-point bending test is expressed follows:  where XE is the distance from the center cross-section of the specimen to an arbitrary element. In Eqs.  and , the stress evaluation point is the center coordinate of each element, and the stress in an element is assumed to be uniform.The following equation was adopted for the fracture mechanics model to evaluate the local fracture stress  around the intrinsic defects:  In the aforementioned equation, F is the geometric factor and c the initial crack length. The most physically reasonable descriptions of the geometric factors in Eq.  for internal defects is expressed as a combination of the stress concentration (pore) and an initial crack[25,50,51]. Based on the fractography, we assumed a fracture mechanics model comprising an ellipsoidal (oblate) pore corresponding to the stress concentration and a circumferential initial crack around the pore, as shown in Figure 2(b)[28–32]. The circumferential initial crack length was assumed to be the length corresponding to one grain around the pore. Hence, the distribution of the grain size c, as presented in Section 3.2, was used.The geometric factor for an initial crack around an ellipsoidal pore, as shown in Figure 2(b), can be expressed as follows[52]:      In the aforementioned equations, Kt is the stress concentration factor and  the notch root radius, which is defined by . Further, . The concrete expression for Kt was presented by Ozaki et al.[28]. We explain how to numerically obtain the bending-strength scatter using the aforementioned two stress evaluation models. Based on the number density of the pores ND, NE number of pores were generated for each element in the specimen model (Figure 2(a)) according to the volume of the element. In this case, the major pore radius R was set randomly based on the GP distribution using an inverse function method based on the cumulative distribution function.Because the largest pore among NE pores in an element can be a candidate for the fracture initiation point of each element, the largest pore RE was extracted for each element as follows:  where  is the randomly generated pore size based on the GP distribution. In addition to RE, the initial crack length c and pore aspect ratio A were set for each element as random numbers following log-normal and normal distributions, respectively, which were applied to the fracture mechanics model of Eqs. –. The aforementioned operations yielded the geometric factor F and local fracture stress  of the potential fracture initiation point for each element.Further, we evaluated the respective bending loads P when the bending stress  was calculated using Eqs.  and , which exceeded the local fracture stress of each element. The lowest bending load was later applied to Eq.  as the peak load of the specimen, and the bending strength  was calculated. In this study, we used fine-grained alumina with limited R-curve behavior, as reported in a previous study[53]. Hence, the bending stress at the initial fracture and bending strength (peak value of the bending stress) obtained from the bending test were assumed to be similar.Based on the aforementioned calculation flow, the bending strengths of  specimens were predicted for the four models presented in Table 1, which were similar to the experiment. The scatter in bending strength was organized using a two-parameter Weibull distribution (Eq. ) and the effective volume (Eq. ). Additionally, the bending strength scatter of the  specimens was simulated to investigate the strength scatter in detail.5. Simulation results and discussion5.1 Specimen size dependencyFigure 11(a) compares the Weibull plots of the bending strengths obtained from the experiments and simulations. The element size and GP distribution character used in the simulation were a cube of 0.25 mm per edge and distribution C (Figures 9(a) and (b)), respectively, thereby resulting in a number of pores per element of . The other parameters used in the simulations are presented in Table 2. The -and m-values obtained from the data in Figure 11(a) are presented in Table 3. The simulation results obtained using distribution C accurately described the experimental Weibull plots under all the test conditions. The m-values in the simulation appear to be slightly larger than those in the experiment. This may be owing to the limited number of test specimens () and the influence of the discretization of the bending stress distribution in the specimen models.Figure 11(b) shows the relationship between the -values and effective volume  for the experimental and simulation results using GP distributions A, B, and C for the major pore radius. The error bars of the -values indicate the upper and lower limits of the bending strength () corresponding to  and 2.3%, respectively. The error bars in  indicate the values estimated from the 90% confidence intervals of the obtained m-values. The size dependence of the -value estimated using Eq.  are represented by the solid and dashed lines, respectively. Figure 11(b) confirms that the -values and overall bending strengths predicted using pore distribution A are higher than those obtained in the experiment. This is because the PWM method focuses on regions with numerous pores, thereby resulting in an underestimation of the probability of the existence of large pores. Conversely, when distribution B was used, the -values and overall bending strengths were lower than those obtained in the experiment. By fitting the higher limit of scatter in the pore distribution measured using X-ray CT, the probability of the existence of large pores was overestimated. In this case, wherein there is some variation in the measured probability density, the probability of existence of large pores lies between distributions A and B. In fact, the numerical results obtained using distribution C can reasonably estimate the -values (Figure 11(b)).5.2 Discretization dependencyExamining the effect of the discretization of the stress evaluation point (element size) is also important for discussing the versatility of the proposed method. The GP distribution C was used for the simulation. To study the discretization dependency in terms of estimating the strength scatter, specimen models were discretized using different sizes of cubic elements with dimensions of 0.125 and 0.5 mm per edge in addition to 0.25 mm. The simulation parameters are presented in Table 2. The number of pores generated per element NE was adjusted based on the element volume. It is noteworthy that the number density of pores ND is identical for all models, irrespective of the discretization size.The Weibull plots created from the simulation results for element sizes of 0.125 and 0.5 mm are shown in Figures 12(a) and (b), respectively. The results for the element size of 0.25 mm are shown in Figure 11(a). Figure 13 shows a summary of the relationship between the strength scatter and effective volume . As shown in Figures 12 and 13, the characteristics of the -values and lower limit of the strength decrease as the effective volume  increases, which are reproduced independently of the element size. Particularly, the results for the element sizes of 0.25 and 0.125 mm were similar, whereas only those of 0.5 mm were slightly higher. The proposed method based on the GP model can obtain similar results even when the element size varies within a certain range.The overall higher predicted strengths for the element size of 0.5 mm are not because of the statistical treatment of the major pore radii using the GP model but because of the accuracy of the calculation of the stress distribution . In the case of simple tensile analysis, even the element size of 0.5 mm yields almost similar results as the element size of 0.125 mm. However, in the bending tests, the finer the discretization, the more accurately the stress distribution can be evaluated (Eqs.  and ). Hence, even when the proposed method is used, sufficient discretization is necessary to reproduce the stress gradient and stress concentration based on the shapes of the components and loading conditions. This restriction is similar to that of ordinary FEA. The proposed method is expected to be applicable for discretizing a single model (specimen or component) by combining elements of different sizes because it uses the number-density information fitted by the GP model.5.3 DiscussionThe proposed numerical simulation method can reasonably estimate the size-dependent bending strength and its scatter directly from the microstructural features, suggesting the various advantages for reliability-assurance of brittle ceramics. We summarized the size dependency of m-values (Figure 14(a)) and -values (Figure 14(b)) as the main indices of scatter and strength, respectively. The error bars of the m-values in Figure 14(a) correspond to 90% confidence intervals[24]. The error bars of the -values indicate the bending strengths corresponding to  and 2.3%. The simulated results obtained from the  bending tests are shown in the figures 14(a) and (b). The error bars of the m-value for  were sufficiently small to be ignored[18,24].As shown in Figure 14(a), the experimental m-values showed a large scatter in the case of  tests, and the m-values were not constant depending on . The simulated m-values also showed some scatter in the  test. Conventionally, a Weibull plot obtained from small and simple-shaped test pieces is used as a reference to estimate the reliability (fracture probability) of ceramic components, assuming that the m-value is constant with respect to the effective volume (Eq. ). Our results imply that it is not possible to determine whether the m-value is independent of  in  tests. More importantly, this indicates that the m-value simulated by the  tests is almost independent of , at least within the tested  range. As shown in Figure 14(b), the size dependence of -values is often estimated using Eq.  within the range of the effective volume, for which the m-value is assumed to be constant. However, the conventional method is required to evaluate the effective volume in advance, ideally in more than 1000 tests, and to obtain information from a reference strength test. In some cases, there is a risk of danger-side reliability evaluations. Contrarily, it is noteworthy that the proposed method does not require an effective volume evaluation in advance and naturally predicts strength scatters and their lower limits for specimens of arbitrary sizes. Another significant advantage of this method is its ability to perform calculations more than 1000 times at high speed and low cost. As shown in Figure 14(b), the proposed method can predict the strength corresponding to an arbitrary fracture probability (99.9, 97.7, 63.2, 10, 2.3, 1, and 0.1%) for an arbitrary number of tests, thereby resulting in a safe-side evaluation of the reliability. Thus, the features of this method provide a significant advantage in ensuring the reliability of brittle ceramics.Generally, ceramic components have larger and more complicated shapes than the test specimens. The proposed method can be applied to the FEA used for the stress analysis of various components with complicated shapes in accordance with previous studies[28–32]. This enables us to examine the design of components in accordance with safety requirements, considering the lower limit of strength even for brittle materials.One open issue is how to fit pore size distribution to the GP model reasonably. In addition to the PWM method, several other methods such as the maximum likelihood method could be used and must be considered in the future. Furthermore, it is necessary to enrich the large-pore data by extending the measurement volume of the X-ray CT and area of SEM. This also makes it possible to extract only large defects composed of agglomerated pores around unsintered powders using a sample mean excess plot.6. ConclusionIn this study, a numerical simulation method for the strength scatter of ceramics based on microstructural data was proposed by employing the GP model, a type of extreme value statistic, to estimate the internal defects that could be the origin of fracture. The validity of the proposed method was verified by comparing its results to the experimental results of four types of bending tests. Based on the results of the X-ray CT measurements, the strength scatter and its lower limit for the same lot of specimens could be predicted by applying the GP model to the probability density function of the major pore radius, regardless of the specimen size and loading condition. Additionally, simulations were conducted using different discretization sizes of the specimens, and the validity of the proposed method for the size dependence of the bending strength and discretization size independence was demonstrated. Unlike conventional metallic materials, it is difficult to evaluate the strength scatter in ceramics based on the results of standardized tests, such as bending tests. However, the proposed method can overcome this problem.Although the R-curve behavior was ignored in the present numerical simulation, it can be handled using the FEA method, which implements a damage model[28–32]. In the future, we plan to expand the proposed method to FEA to consider different microstructures, strengths, toughness, and R-curve behaviors.AcknowledgementsWe express our gratitude to Dr. Makoto Watanabe and Ms. Akiko Takenouchi for their experimental assistance and technical guidance regarding the X-ray CT facilities. We also thank Ms. Mariko Iguchi for her assistance with the electron microscopic analysis. Finally, we thank Editage (www.editage.jp) for the English language editing. Funding:This study is based on the results obtained from project JPNP22005, commissioned by the New Energy and Industrial Technology Development Organization (NEDO). AppendixIn this Appendix, we describe in detail the fitting of the GP model to the pore size distribution of AS999.Customarily, we determined the threshold value as the minimum value when plots to the right of the threshold value were considered close to a straight line. Hence, an optimal threshold must be selected for accurate data analysis. When the threshold is small, more data points exceed it, thereby worsening the fit of the GP distribution. Hence, the variance in the estimation reduced, although the bias increased. Conversely, when the threshold was large, there were fewer data points. However, the GP distribution was better fitted. In this case, the estimation bias decreased, whereas the variance increased. As shown in Figure 8, the threshold was set at 5.9 μm, while the data to the right of the threshold were considered to contain the straight line with the positive slope. In the data with 1399 pore sizes, 891 data points were above the threshold value of 5.9 μm. It is noteworthy that the number of data above a larger threshold (14 μm), which appears to be a straight line at first impression, is only 32.Figures 15(a) and (b) show the probability and quantile plots, respectively, based on the GP distribution, wherein the threshold value was 5.9 μm. The probability plot shows that the data were distributed without a considerable offset from the diagonal line, thereby indicating that the fit of the GP distribution was sufficient for the pore size distribution in the target ceramics. The quantile plot also showed that the data were distributed without deviating from a linear relationship, thereby confirming the fit to the GP distribution. Only the rightmost data deviated from the straight line. However, it was within the 95% confidence interval and considered to be within the range of variation caused by the small amount of data.References[1] Y. Sakabe, Multilayer ceramic capacitors, Curr Opin Solid State Mater Sci. 2 (1997) 584–587.[2] P. Zhao, Z. Cai, L. Wu, C. Zhu, L. Li, X. Wang, Perspectives and challenges for lead-free energy-storage multilayer ceramic capacitors, Journal of Advanced Ceramics. 10 (2021) 1153–1193. https://doi.org/10.1007/s40145-021-0516-8.[3] Z. Lu, G. Wang, W. Bao, J. Li, L. Li, A. Mostaed, H. Yang, H. Ji, D. Li, A. Feteira, F. Xu, D.C. Sinclair, D. Wang, S.Y. Liu, I.M. Reaney, Superior energy density through tailored dopant strategies in multilayer ceramic capacitors, Energy Environ Sci. 13 (2020) 2938–2948. https://doi.org/10.1039/d0ee02104k.[4] H. Zhang, T. Wei, Q. Zhang, W. Ma, P. Fan, D. Salamon, S.T. Zhang, B. Nan, H. Tan, Z.G. Ye, A review on the development of lead-free ferroelectric energy-storage ceramics and multilayer capacitors, J Mater Chem C Mater. 8 (2020) 16648–16667. https://doi.org/10.1039/d0tc04381h.[5] J. Robert Kelly, Dental ceramics: Current thinking and trends, Dental Clinics. 48 (2004) 513–530. https://doi.org/10.1016/j.cden.2004.01.003.[6] I. Denry, J.A. Holloway, Ceramics for dental applications: A review, Materials. 3 (2010) 351–368. https://doi.org/10.3390/ma3010351.[7] J. Han, F. Zhang, B. Van Meerbeek, J. Vleugels, A. Braem, S. Castagne, Laser surface texturing of zirconia-based ceramics for dental applications: A review, Materials Science and Engineering C. 123 (2021) 112034. https://doi.org/10.1016/j.msec.2021.112034.[8] A.N. Samant, N.B. Dahotre, Laser machining of structural ceramics-A review, J Eur Ceram Soc. 29 (2009) 969–993. https://doi.org/10.1016/j.jeurceramsoc.2008.11.010.[9] K.N. Lee, Current status of environmental barrier coatings for Si-Based ceramics, Surf Coat Technol. 133 (2000) 1–7.[10] H. Ohnabe, S. Masaki, M. Onozuka, K. Miyahara, T. Sasa, Potential application of ceramic matrix composites to aero-engine components, Compos Part A Appl Sci Manuf. 30 (1999) 489–496.[11] O. Gavalda Diaz, G. Garcia Luna, Z. Liao, D. Axinte, The new challenges of machining Ceramic Matrix Composites (CMCs): Review of surface integrity, Int J Mach Tools Manuf. 139 (2019) 24–36. https://doi.org/10.1016/j.ijmachtools.2019.01.003.[12] A.G. EVANS, Structural Reliability:A Processing‐Dependent Phenomenon, Journal of the American Ceramic Society. 65 (1982) 127–137. https://doi.org/10.1111/j.1151-2916.1982.tb10380.x.[13] H. Peterlik, Relationship of Strength and Defects of Ceramic Materials and Their Treatment by Weibull Theory, Journal of the Ceramic Society of Japan. 109 (2001) S121–S126.[14] A. Wang, P. Hu, X. Zhao, Z. Wang, C. Zhang, Y. Wang, Modelling and experimental investigation of pore-like flaw-strength response in structural ceramics, Ceram Int. 46 (2020) 14431–14438. https://doi.org/10.1016/j.ceramint.2020.02.240.[15] G.D. Quinn, R. Morrell, Design Data for Engineering Ceramics: A Review of the Flexure Test, Journal of the American Ceramic Society. 74 (1991) 2037–2066. https://doi.org/10.1111/j.1151-2916.1991.tb08259.x.[16] R. Danzer, A General Strength Distribution Function for Brittle Materials, J Eur Ceram Soc. 10 (1992) 461–472.[17] R. Danzer, P. Supancic, J. Pascual, T. Lube, Fracture statistics of ceramics - Weibull statistics and deviations from Weibull statistics, Eng Fract Mech. 74 (2007) 2919–2932. https://doi.org/10.1016/j.engfracmech.2006.05.028.[18] R. Danzer, T. Lube, P. Supancic, R. Damani, Fracture of ceramics, Adv Eng Mater. 10 (2008) 275–298. https://doi.org/10.1002/adem.200700347.[19] J.B. Quinn, G.D. Quinn, A practical and systematic review of Weibull statistics for reporting strengths of dental materials, Dental Materials. 26 (2010) 135–147. https://doi.org/10.1016/j.dental.2009.09.006.[20] B. Deng, D. Jiang, J. Gong, Is a three-parameter Weibull function really necessary for the characterization of the statistical variation of the strength of brittle ceramics?, J Eur Ceram Soc. 38 (2018) 2234–2242. https://doi.org/10.1016/j.jeurceramsoc.2017.10.017.[21] Y. Lu, Z. Mei, J. Zhang, S. Gao, X. Yang, B. Dong, L. Yue, H. Yu, Flexural strength and Weibull analysis of Y-TZP fabricated by stereolithographic additive manufacturing and subtractive manufacturing, J Eur Ceram Soc. 40 (2020) 826–834. https://doi.org/10.1016/j.jeurceramsoc.2019.10.058.[22] S. Nohut, Three-parameter (3P) weibull distribution for characterization of strength of ceramics showing R-Curve behavior, Ceram Int. 47 (2021) 2270–2279. https://doi.org/10.1016/j.ceramint.2020.09.067.[23] Z.P. Bažant, Size effect on structural strength: a review, Archive of Applied Mechanics. 69 (1999) 703–725.[24] ASTM C1239-07, ASTM Annual Book of Standards, vol 15.01. American Society for Testing and Materials, West Conshohocken. (2009). https://doi:10.1520/C1239-07.[25] A. Zimmermann, J. Rödel, Fracture statistics based on pore/grain-size interaction, Journal of the American Ceramic Society. 82 (1999) 2279–2281. https://doi.org/10.1111/j.1151-2916.1999.tb02080.x.[26] S. Nakamura, S. Tanaka, Z. Kato, K. Uematsu, Strength-processing defects relationship based on micrographic analysis and fracture mechanics in alumina ceramics, Journal of the American Ceramic Society. 92 (2009) 688–693. https://doi.org/10.1111/j.1551-2916.2008.02904.x.[27] S. Nohut, Influence of sample size on strength distribution of advanced ceramics, Ceram Int. 40 (2014) 4285–4295. https://doi.org/10.1016/j.ceramint.2013.08.093.[28] S. Ozaki, Y. Aoki, T. Osada, K. Takeo, W. Nakao, Finite element analysis of fracture statistics of ceramics: Effects of grain size and pore size distributions, Journal of the American Ceramic Society. 101 (2018) 3191–3204. https://doi.org/10.1111/jace.15468.[29] S. Ozaki, M. Nakamura, T. Osada, Finite element analysis of the fracture statistics of self-healing ceramics, Sci Technol Adv Mater. 21 (2020) 609–625. https://doi.org/10.1080/14686996.2020.1800368.[30] C. Ito, T. Osada, S. Ozaki, Finite element analysis of fracture behavior in ceramics: Competition between artificial notch and internal defects under three-point bending, Ceram Int. 48 (2022) 36460–36468. https://doi.org/10.1016/j.ceramint.2022.08.206.[31] K. Takeo, Y. Aoki, T. Osada, W. Nakao, S. Ozaki, Finite element analysis of the size effect on ceramic strength, Materials. 12 (2019) 2885. https://doi.org/10.3390/ma12182885.[32] S. Ozaki, K. Yamagata, C. Ito, T. Kohata, T. Osada, Finite element analysis of fracture behavior in ceramics: Prediction of strength distribution using microstructural features, Journal of the American Ceramic Society. 105 (2022) 2182–2195. https://doi.org/10.1111/jace.18201.[33] G. Okuma, T. Osada, H. Minagawa, Y. Arai, R. Inoue, H. Kakisawa, K. Shimoda, A. Takeuchi, M. Uesugi, S. Tanaka, F. Wakai, Heterogeneities and defects in powder compacts and sintered alumina bodies visualized by using the synchrotron X-ray CT, J Eur Ceram Soc. 43 (2023) 486–492. https://doi.org/10.1016/j.jeurceramsoc.2022.10.020.[34] G. Okuma, M. Endo, H. Minagawa, R. Inoue, H. Kakisawa, T. Kohata, T. Osada, T. Yamamoto, M. Azuma, A. Takeuchi, M. Uesugi, O. Guillon, F. Wakai, 3D Visualization of Morphological Evolution of Large Defects during Spark Plasma Sintering of Alumina Granules, Adv Eng Mater. (2023) 2201534. https://doi.org/10.1002/adem.202201534.[35] A. Wormsen, B. Sjödin, G. Härkegård, A. Fjeldstad, Non-local stress approach for fatigue assessment based on weakest-link theory and statistics of extremes, Fatigue Fract Eng Mater Struct. 30 (2007) 1214–1227. https://doi.org/10.1111/j.1460-2695.2007.01190.x.[36] A. Fjeldstad, A. Wormsen, G. Härkegård, Simulation of fatigue crack growth in components with random defects, Eng Fract Mech. 75 (2008) 1184–1203. https://doi.org/10.1016/j.engfracmech.2007.04.006.[37] A. Wormsen, A. Fjeldstad, G. Härkegård, A post-processor for fatigue crack growth analysis based on a finite element stress field, Comput Methods Appl Mech Eng. 197 (2008) 834–845. https://doi.org/10.1016/j.cma.2007.09.012.[38] S. Beretta, More than 25 years of extreme value statistics for defects: Fundamentals, historical developments, recent applications, Int J Fatigue. 151 (2021) 106407. https://doi.org/10.1016/j.ijfatigue.2021.106407.[39] S. Beretta, C. Anderson, Y. Murakami, Extreme value models for the assessment of steels containing multiple types of inclusion, Acta Mater. 54 (2006) 2277–2289. https://doi.org/10.1016/j.actamat.2006.01.016.[40] L. ‐Y Chao, D.K. Shetty, Extreme‐Value Statistics Analysis of Fracture Strengths of a Sintered Silicon Nitride Failing from Pores, Journal of the American Ceramic Society. 75 (1992) 2116–2124. https://doi.org/10.1111/j.1151-2916.1992.tb04473.x.[41] Z.P. Bažant, S.D. Pang, Activation energy based extreme value statistics and size effect in brittle and quasibrittle fracture, J Mech Phys Solids. 55 (2007) 91–131. https://doi.org/10.1016/j.jmps.2006.05.007.[42] J. Pickands III, Statistical Inference Using Extreme Order Statistics, The Annals of Statistics. 3 (1975) 119–131.[43] N.A. WEIL, I.M. DANIEL, Analysis of fracture probabilities in nonuniformly stressed brittle materials, Journal of the American Ceramic Society. 47 (1964) 268–274. https://doi.org/10.1111/j.1151-2916.1964.tb14413.x.[44] K. Ando, Y. Shirai, M. Nakatani, Y. Kobayashi, S. Sato, (Crack-healing+proof test): a new methodology to guarantee the structural integrity of a ceramics component, J Eur Ceram Soc. 22 (2002) 121–128. www.elsevier.com/locate/jeurceramsoc.[45] T. Osada, T. Hara, M. Mitome, S. Ozaki, T. Abe, K. Kamoda, T. Ohmura, Self-healing by design: universal kinetic model of strength recovery in self-healing ceramics, Sci Technol Adv Mater. 21 (2020) 593–608. https://doi.org/10.1080/14686996.2020.1796468.[46] J. Takahashi, H. Suito, Evaluation of the Accuracy of the Three-Dimensional Size Distribution Estimated from the Schwartz-Saltykov Method, Metallurgical and Materials Transactions A. 34 (2003) 171–181.[47] S.A. Saltykov, Stereometric Metallography (2nd edn.), Metallurgizdat, Moscow, 1958.[48] T. ‐S Yeh, M.D. Sacks, Effect of Particle Size Distribution on the Sintering of Alumina, Journal of the American Ceramic Society. 71 (1988) C‐484-C‐487. https://doi.org/10.1111/j.1151-2916.1988.tb05812.x.[49] J.R.M. Hosking, J.R. Wallis, Parameter and Quantile Estimation for the Generalized Pareto Distribution, Technometrics. 29 (1987) 339–349.[50] D.J. Green, Stress intensity factor estimates for annular cracks at spherical voids, Journal of the American Ceramic Society. 63 (1980) 342–344. https://doi.org/10.1111/j.1151-2916.1980.tb10737.x.[51] A. Zimmermann, J. Rödel, Generalized Orowan‐Petch Plot for Brittle Fracture, Journal of the American Ceramic Society. 81 (1998) 2527–2532.[52] Y. Murakami, Stress Intensity Factor Handbook 1(5.4), Oxford: Pergamon Press. (1987) 249–254.[53] H. Tomaszewski, M. Boniecki, H. Weglarz, Effect of grain size on R-curve behaviour of alumina ceramics, J Eur Ceram Soc. 20 (2000) 2569–2574.Figure and table captionsFig. 1.  Application of extreme value statistics to the probability density function of large pores, which could be the fracture initiation points: (a) schematic of the size distributions of three types of pores and the extraction of pore distribution of possible fracture initiation points using the Generalized Pareto model; (b) schematic of the different types of pores with different generation mechanisms.Fig. 2.  Bending test analysis model with the fracture mechanics model implemented: (a) discretization of the specimen using cubic elements and the microstructural data set for each element; (b) schematic of the fracture mechanics model (the circumferential circular crack emanating from an ellipsoidal pore).Fig. 3.  2D images obtained using SEM for alumina ceramics (AS999): (a) bending test specimen geometry and observed area on the cross section of the specimen, (b) number of SEM images on the cross section, and (c) high-resolution image showing the types of internal pores classified using the machine-learning function in the ImageJ software.Fig. 4.  SEM image of the fracture surface of a specimen after the bending test (3pb-S, bending strength: 335 MPa). Here, a large pore including unsintered powders is confirmed at the fracture initiation point.Fig. 5.  3D image obtained using X-ray CT for alumina ceramics (AS999). The image shows 1399 pores within a volume of 220.52×214.02×279.85 µm.Fig. 6.  Size distributions of different types of pores obtained from SEM images: (a) 2D size distribution of each type of pore (small pores and large pores) and their combination; (b) comparison of pore size distribution per unit volume converted using the SS method.Fig. 7.  Grain size features: (a) EBSD images; (b) cumulative frequency of the grain size distribution measured using EBSD and the best-fit curves using the log-normal distribution with an expectation and variance of Ec = 6.8×10-3 mm and Vc = 12.6×10-6 mm2, respectively. Fig. 8.  Sample mean excess plot of pore size obtained using X-ray CT. Fig. 9.  Approximation of the pore size distribution using the GP distribution: (a) probability density functions; (b) exponential probability plots.Fig. 10  Cumulative frequency of the aspect ratio obtained from pores with sizes larger than 5.9 μm and the best-fit curves using the normal distribution with a mean value and variance of A = 0.485 and σA = 0.156, respectively.Fig. 11. Comparison of experimental and simulation results of bending strength: (a) comparison of Weibull plots using distribution C; (b) specimen size dependence of bending strength. Three kinds of GP distributions of the pore size are examined, as shown in Figure 9. Closed and open plots correspond to β-values obtained from the experiment and simulation, respectively, and the respective error bars show fracture probabilities of 2.3% and 97.7%. The error bars of  indicate the values estimated from 90% confidence intervals of the obtained m-values. For discrimination purposes, plot sizes are increased in the order of 3pb-S, 3pb-L, 4pb-S, and 4pb-L.Fig. 12. Comparison of Weibull plots between simulation results under different discretization sizes and experimental results. In this case, distribution C is adopted: (a) the discretization size is 0.125 mm; (b) the discretization size is 0.5 mm. Fig. 13. Comparison of experimental and simulation results for the specimen size dependence of bending strength. Three kinds of discretization sizes are examined using distribution C. Closed and open plots correspond to β-values obtained from the experiment and simulation, respectively, and the respective error bars show fracture probabilities of 2.3% and 97.7%. The error bars of  indicate the values estimated from 90% confidence intervals of the obtained m-values. For discrimination purposes, plot sizes are increased in the order of 3pb-S, 3pb-L, 4pb-S, and 4pb-L.Fig. 14. Simulation results of parameters of Weibull plot obtained from N = 1000 specimens using distribution C. Closed and open plots correspond to parameters obtained from the experiment and simulation, respectively: (a) comparison of m-values with experimental results; the error bars of  indicate the values estimated from 90% confidence intervals of the experimentally obtained m-values, and (b) specimen size dependence of strength scatter. The error bars of the experimental data show fracture probabilities of 2.3% and 97.7%. For discrimination purposes, plot sizes are increased in the order of 3pb-S, 3pb-L, 4pb-S, and 4pb-L.Fig. 15. Fitness of the GP model to the pore size distribution: (a) probability plot; (b) quantile plot.Table 1. Dimension of the test specimen and span length of the bending test.Table 2. Parameters of pore distribution used to examine the effect of the discretization size. The number density of pores, ND, is 69151 mm-3.Table 3. - and m-values obtained from the experiments and simulations (N = 30). 40image2.wmf12B23(),2PSStWs-=image44.wmfboleObject50.binimage45.wmf30N=oleObject51.binimage46.wmfboleObject52.binimage47.wmfeffVoleObject53.binoleObject54.binoleObject55.binoleObject2.binimage48.wmfB()97.7%Fs=oleObject56.binoleObject57.binoleObject58.binoleObject59.binoleObject60.binoleObject61.binimage49.wmfeffVoleObject62.binimage50.wmfbimage3.wmfBB()1exp,mFssbìüæö=--íýç÷èøîþoleObject63.binoleObject64.binimage51.wmfEsoleObject65.binoleObject66.binoleObject67.binoleObject68.binimage52.wmf1000N=oleObject69.binoleObject70.binoleObject3.binimage53.wmf30N=oleObject71.binimage54.wmfeffVoleObject72.binoleObject73.binoleObject74.binoleObject75.binimage55.wmf1000N=oleObject76.binoleObject77.binimage4.wmfboleObject78.binoleObject79.binoleObject4.binoleObject5.binimage5.wmfB()FsoleObject6.binimage56.wmfeffVoleObject80.binoleObject81.binoleObject82.binimage6.wmf30N=oleObject7.binimage7.wmfeffVoleObject8.binimage8.wmf212effs11112(1)1,SSSVVmSmSéù-=+êú++ëûoleObject9.binimage9.wmfs1VStW=oleObject10.binoleObject11.binoleObject12.binoleObject13.binimage10.wmfeff00,mmVVbb=oleObject14.binimage11.wmf0boleObject15.binoleObject16.binoleObject17.binimage12.wmf4.0MPamoleObject18.binimage13.wmf11()(,)(),kVAiNjijNia==DåoleObject19.binimage14.wmfmax/DkD=oleObject20.binimage15.wmfmaxDoleObject21.binimage16.wmf(,)ijaoleObject22.binimage17.wmf22221(,)(1).ijjijia=----oleObject23.binimage18.wmf30k=oleObject24.binimage19.wmf1/11for0,()01expfor.RuHRuRuxxxsxs-ì-æö-+ï¹ç÷ïèø-=í-æöï=--ç÷ïèøîoleObject25.binimage20.wmf1/110for1,()10forexp.RuhRuRuxxxssxss--ì-æö¹+ïç÷ïèø-=í-æöï=-ç÷ïèøîoleObject26.binimage21.wmfsoleObject27.binimage22.wmfxoleObject28.binimage23.wmfmax11,():,uniiuuRuuRn=ìüæöïï-<íýç÷ïïèøîþåoleObject29.binimage24.wmf3D69151mmN-=oleObject30.binimage25.wmfEsoleObject31.binimage26.wmfE12ER(),4YSSPIs-=oleObject32.binimage27.wmfE1EER(/2),2YSXPIs-=oleObject33.binimage28.wmfFsoleObject34.binimage29.wmfICF.KFcsp=oleObject35.binimage30.wmf112for/1,max(,)for/1.FcFFFcrr£ì=í>îoleObject36.binimage31.wmf21t2411131.125(10.22380.1643),36(1)(1)FKlllléùìü=+++-êúíý++êúîþëûoleObject37.binimage32.wmf22.RcFcp+=oleObject38.binimage33.wmfroleObject39.binimage1.wmfBsimage34.wmf2b/RRr=oleObject40.binimage35.wmf/clr=oleObject41.binimage36.wmfEGPDGPDGPDGPDE123max(,,,,),NRRRRR=LoleObject42.binimage37.wmfGPDiRoleObject43.binimage38.wmfFsoleObject44.binoleObject1.binimage39.wmfEsoleObject45.binimage40.wmfBsoleObject46.binimage41.wmf30N=oleObject47.binimage42.wmf1000N=oleObject48.binimage43.wmfE1080N=oleObject49.bin