Probability Density Function Models for Float Glass under Mechanical Loading with Varying Parameters

Glass as a construction material has become indispensable and is still on the rise in the building industry. However, there is still a need for numerical models that can predict the strength of structural glass in different configurations. The complexity lies in the failure of glass elements largely driven by pre-existing microscopic surface flaws. These flaws are present over the entire glass surface, and the properties of each flaw vary. Therefore, the fracture strength of glass is described by a probability function and will depend on the size of the panels, the loading conditions and the flaw size distribution. This paper extends the strength prediction model of Osnes et al. with the model selection by the Akaike information criterion. This allows us to determine the most appropriate probability density function describing the glass panel strength. The analyses indicate that the most appropriate model is mainly affected by the number of flaws subjected to the maximum tensile stresses. When many flaws are loaded, the strength is better described by a normal or Weibull distribution. When few flaws are loaded, the distribution tends more towards a Gumbel distribution. A parameter study is performed to examine the most important and influencing parameters in the strength prediction model.


Introduction
In the building industry, soda-lime-silica glass has been mostly used as infill panels, but now more and more commonly as a structural component. As such, one of its most significant challenges is the brittle behavior, resulting in an instantaneous fracture. This challenge can be partly tackled by lamination or combining the glass with, e.g., steel reinforcement [1][2][3][4][5]. However, predicting the behavior of a structural glass panel and its corresponding failure strength is not easy because of the presence and randomness of pre-existing surface flaws on the glass surface, known as 'Griffith flaws' [6][7][8][9]. These flaws initiate the failure of the glass once its maximum tensile strength is exceeded. They are created during production, transportation, cutting, finishing, etc., and their properties vary from flaw to flaw. Consequently, the fracture strength of the glass panels may be described by a probability function, which will depend on the size of the glass panels and the loading conditions. Due to this complexity, glass plates are frequently modeled with a deterministic fracture strength, making the calculation easier. However, this approach does not properly predict the failure strength of glass and does not allow for an explicit evaluation of the safety level (e.g., to allow for a reliability-based design format). Therefore, models need to predict the probabilistic failure strength of the glass in all kinds of configurations.
Some models are already available in the literature, going from very advanced models to more simplified models [10][11][12]. One of the latter is a strength prediction model of annealed glass proposed by Yankelevsky [12] and later adopted by Osnes [13]. The analysis obtains the strength of glass panels as a result and not as an assumption, which is an advantage compared to advanced models. This model predicts the failure of glass elements Table 1. Some thermal and mechanical characteristics of annealed soda-lime-silica glass [17][18][19].

Strength Prediction Model
A crack in a glass pane is initiated when the Mode I stress intensity factor K I of a flaw reaches the fracture toughness K Ic [16]. This crack initiation criterion is presented in Equation (1). The stress intensity factor is calculated with Equation (2), where Y stands for the geometry factor, σ is the stress perpendicular to the flaw orientation and a represents the flaw depth [16]. The geometry factor Y is a dimensionless factor that is dependent on the geometry of the flaw. Further background on the calculation of Y and the other parameters in this crack initiation criterion is given in Section 2.2.2. Glass panes will only fail if flaws on their surface are subjected to tension perpendicular to the flaw orientation. To numerically simulate the crack initiation, surface flaws must first be distributed over the glass surfaces, each with a flaw orientation α i and a flaw depth a i . Once these flaw orientations and depths are generated, they are distributed over the glass surfaces. Next, a mechanical load is applied, which induces stresses. From these stresses, the normal stress σ for each flaw (i.e., the stresses perpendicular to the flaw orientation) can be calculated by Equation (3). The stresses σ x , σ y and τ xy can be calculated using a finite element (FE) program such as Abaqus [20] or calculated analytically for simple configurations.
Once the normal stresses are calculated, Equation (2) can be evaluated at a specific flaw location considering the corresponding flaw depth a i and a geometry factor Y. This results in a stress intensity factor K I for each flaw. Subsequently, the crack initiation criterion of Equation (1) is evaluated for each flaw. If the corresponding stress intensity factor K I exceeds the fracture toughness K Ic , the fracture is initiated in that particular flaw, resulting in fracture of the total glass panel. On the other hand, if the crack initiation criterion is not met, the mechanical load is increased, resulting in increased stresses σ and increased stress intensity factors K I to eventually check again if one of these stress intensity factors K I exceeds the fracture toughness K Ic . This procedure is repeated until the glass eventually fractures, resulting in an evaluation of the maximum load for the glass pane. A flowchart of this procedure is presented in Figure 1.
the other parameters in this crack initiation criterion is given in Section 2.2.2.

Procedure
Glass panes will only fail if flaws on their surface are subjected to tension perpendicular to the flaw orientation. To numerically simulate the crack initiation, surface flaws must first be distributed over the glass surfaces, each with a flaw orientation and a flaw depth . Once these flaw orientations and depths are generated, they are distributed over the glass surfaces. Next, a mechanical load is applied, which induces stresses. From these stresses, the normal stress for each flaw (i.e., the stresses perpendicular to the flaw orientation) can be calculated by Equation (3). The stresses , and can be calculated using a finite element (FE) program such as Abaqus [20] or calculated analytically for simple configurations.
Once the normal stresses are calculated, Equation (2) can be evaluated at a specific flaw location considering the corresponding flaw depth and a geometry factor . This results in a stress intensity factor for each flaw. Subsequently, the crack initiation criterion of Equation (1) is evaluated for each flaw. If the corresponding stress intensity factor exceeds the fracture toughness , the fracture is initiated in that particular flaw, resulting in fracture of the total glass panel. On the other hand, if the crack initiation criterion is not met, the mechanical load is increased, resulting in increased stresses and increased stress intensity factors to eventually check again if one of these stress intensity factors exceeds the fracture toughness . This procedure is repeated until the glass eventually fractures, resulting in an evaluation of the maximum load for the glass pane. A flowchart of this procedure is presented in Figure 1.

Parameters
The SPM relies on the dimensions of the setup and six other parameters: the fracture toughness, the flaw density, the flaw orientation, the flaw depth, the maximum flaw depth, and the flaw shape. The fracture toughness, at inert conditions, K Ic is fixed for an elliptical flaw as 0.75 MPa √ m [18,21].
To quantify the number of flaws simulated on the glass surfaces, the flaw density parameter ρ f law is used, representing the number of flaws present in one square centimeter. The literature shows that two flaws for each square centimeter are a good assumption [13,22].
With this information, the next two parameters can be calculated: the flaw orientation α i and the flaw depth a i for each flaw. The former represents a random value uniformly distributed between zero and π, and the latter represents a random value between zero and the maximum flaw depth a max .
It may be assumed that the flaws distribution function is of an exponential shape which can be expressed with the following Equation (4) [12,23]: where N 0 represents the number of surface flaws present on a jumbo glass plate (3210 mm by 6000 mm), N i the number of flaws that are larger or equal to the maximum flaw depth a max , a i the size of a certain flaw on the glass plate, and η the distribution parameter. Suppose a i = a max ; it means N i will be equal to 1. Substituting this condition in Equation (4), one obtains the distribution parameter η equal to a max /ln(N 0 ) [12]. With this information, Equation (4) can be rewritten to determine the random depth value a i for each flaw as the following Equation (5): where N i is selected with the expression N i = R(N 0 − 1) + 1 in which R is a random value uniformly distributed between zero and one. The maximum flaw depth parameter a max is considered equal to 0.1 mm [13].
Next, information about the geometry of the flaw is needed to calculate the geometry factor Y. However, these flaws can have various shapes, so some assumptions must be made until more information is available. Therefore, it is chosen to model the flaws on the glass surface with an elliptical shape [12], as depicted in Figure 2. The model assumes the flaws are located in the center of each square centimeter of the surfaces. Levengood [23] stated that for Griffith flaws, one-half of the crack length c was observed to be of the same order of magnitude as the crack depth a. Based on this, the parameter of flaw depth on half length of the flaw a/c for elliptical flaws is considered 1 and is called the flaw shape.

Parameters
The SPM relies on the dimensions of the setup and six other parameters: the fracture toughness, the flaw density, the flaw orientation, the flaw depth, the maximum flaw depth, and the flaw shape. The fracture toughness, at inert conditions, is fixed for an elliptical flaw as 0.75 MPa√m [18,21].
To quantify the number of flaws simulated on the glass surfaces, the flaw density parameter is used, representing the number of flaws present in one square centimeter. The literature shows that two flaws for each square centimeter are a good assumption [13,22].
With this information, the next two parameters can be calculated: the flaw orientation and the flaw depth for each flaw. The former represents a random value uniformly distributed between zero and π, and the latter represents a random value between zero and the maximum flaw depth . It may be assumed that the flaws distribution function is of an exponential shape which can be expressed with the following Equation (4) [12,23]: where 0 represents the number of surface flaws present on a jumbo glass plate (3210 mm by 6000 mm), the number of flaws that are larger or equal to the maximum flaw depth , the size of a certain flaw on the glass plate, and the distribution parameter. Suppose = ; it means will be equal to 1. Substituting this condition in Equation (4), one obtains the distribution parameter equal to /ln( 0 ) [12]. With this information, Equation (4) can be rewritten to determine the random depth value for each flaw as the following Equation (5): where is selected with the expression = ( − 1) + 1 in which is a random value uniformly distributed between zero and one. The maximum flaw depth parameter is considered equal to 0.1 mm [13].
Next, information about the geometry of the flaw is needed to calculate the geometry factor . However, these flaws can have various shapes, so some assumptions must be made until more information is available. Therefore, it is chosen to model the flaws on the glass surface with an elliptical shape [12], as depicted in Figure 2. The model assumes the flaws are located in the center of each square centimeter of the surfaces. Levengood [23] stated that for Griffith flaws, one-half of the crack length was observed to be of the same order of magnitude as the crack depth . Based on this, the parameter of flaw depth on half length of the flaw / for elliptical flaws is considered 1 and is called the flaw shape.  The geometry factor Y of the elliptical flaws considered here is calculated using Equation (6), where λ s is the surface correction factor, f (∅) is the angular function and Q is the flaw shape parameter [6,13,16]. All the parameters in Equation (6) are calculated using the empirical expressions (7) to (9), with ∅ equal to zero, which stands for the angle of the point on the elliptical flaw at which the stress intensity factor is at its maximum. Taking this into account, the final geometry factor Y for elliptical flaws, calculated with Equation (6), is equal to 0.729 [6,13,16].

Case Study
To demonstrate the different steps in the SPM, a case study performed by Osnes et al. [13] will be used. The experimental data of this case study will be used to validate the outcome of the SPM at the end of this section.
The case study is a four-point bending test performed on glass specimens of 300 mm long, 60 mm wide and 4 mm thick, where the support span L s is equal to 280 mm and the loading span L l is equal to half of the support span, i.e., 140 mm. The dimensions of the setup and the specimens are based on the ASTM standard C1161-13 [24], originally intended for testing advanced ceramics.

SPM Applied on Case Study
The above-mentioned four-point bending setup is modeled in Python. The area of the glass panel combined with the flaw density ρ f law results in the total amount of flaws present on one glass surface. With this information, the flaws can be randomly distributed over both the top and bottom surfaces. Each flaw can be assigned a random flaw orientation α i and a random flaw depth a i as discussed in Section 2.2.2. To demonstrate this concept, Figures 3 and 4 visualize the distribution of the flaw orientations and depths, respectively. To display this clearly, these visualizations are divided into elements of one square centimeter and exceptionally adopt the simplification of displaying one color for each square centimeter, i.e., a flaw density of one flaw/cm 2 .
The geometry factor of the elliptical flaws considered here is calculated using Equation (6), where is the surface correction factor, (∅) is the angular function and is the flaw shape parameter [6,13,16].
All the parameters in Equation (6) are calculated using the empirical expressions (7) to (9), with ∅ equal to zero, which stands for the angle of the point on the elliptical flaw at which the stress intensity factor is at its maximum. Taking this into account, the final geometry factor for elliptical flaws, calculated with Equation (6), is equal to 0.729 [6,13,16].

Case Study
To demonstrate the different steps in the SPM, a case study performed by Osnes et al. [13] will be used. The experimental data of this case study will be used to validate the outcome of the SPM at the end of this section.
The case study is a four-point bending test performed on glass specimens of 300 mm long, 60 mm wide and 4 mm thick, where the support span is equal to 280 mm and the loading span is equal to half of the support span, i.e., 140 mm. The dimensions of the setup and the specimens are based on the ASTM standard C1161-13 [24], originally intended for testing advanced ceramics.

SPM Applied on Case Study
The above-mentioned four-point bending setup is modeled in Python. The area of the glass panel combined with the flaw density results in the total amount of flaws present on one glass surface. With this information, the flaws can be randomly distributed over both the top and bottom surfaces. Each flaw can be assigned a random flaw orientation and a random flaw depth as discussed in Section 2.2.2. To demonstrate this concept, Figures 3 and 4 visualize the distribution of the flaw orientations and depths, respectively. To display this clearly, these visualizations are divided into elements of one square centimeter and exceptionally adopt the simplification of displaying one color for each square centimeter, i.e., a flaw density of one flaw/cm 2 .   Next, the mechanical load is applied which will induce stresses. As mentioned, these stresses can be calculated using an FE program such as Abaqus or calculated analytically. For visualizing purposes, Abaqus is used here to calculate the stresses.
Both the four-point bending setup and the increased mechanical loading are numerically simulated. This simulation's output are the stress states induced by the mechanical loading at evenly spaced loading intervals and for every location on the glass pane. Stress states consist of the in-plane normal stresses in the x and y direction, respectively, called σ x and σ y , and the in-plane shear stress τ xy . In Figure 5, the σ x stress states of the bottom and top surfaces are visualized for a total load of 1000 N and 1800 N.  Next, the mechanical load is applied which will induce stresses. As mentioned, these stresses can be calculated using an FE program such as Abaqus or calculated analytically. For visualizing purposes, Abaqus is used here to calculate the stresses.
Both the four-point bending setup and the increased mechanical loading are numerically simulated. This simulation's output are the stress states induced by the mechanical loading at evenly spaced loading intervals and for every location on the glass pane. Stress states consist of the in-plane normal stresses in the x and y direction, respectively, called and , and the in-plane shear stress . In Figure 5, the stress states of the bottom and top surfaces are visualized for a total load of 1000 N and 1800 N. These , and stress states are used to calculate the normal stress of each flaw for any load on the structure, as demonstrated in Equation (3). Combining this stress evolution with the corresponding flaw depths and the geometry factor as calculated in Equation (6), the evolution of the stress intensity factor for every flaw can be calculated and visualized, as depicted in Figure 6. The evolution of the stress intensity factor will show the location where a first flaw reaches a stress intensity factor equal to the fracture toughness of 0.75 MPa√m. Once this failing flaw is found, one knows the fracture load and the fracture location for this particular glass pane in a four-point bending setup. Figure 7 visualizes the fracture location for the random glass pane shown in the former figures. In this example, the fracture originated at the bottom surface location with and coordinates equal to 175 mm and 35 mm, respectively.  Next, the mechanical load is applied which will induce stresses. As mentioned, these stresses can be calculated using an FE program such as Abaqus or calculated analytically. For visualizing purposes, Abaqus is used here to calculate the stresses.
Both the four-point bending setup and the increased mechanical loading are numerically simulated. This simulation's output are the stress states induced by the mechanical loading at evenly spaced loading intervals and for every location on the glass pane. Stress states consist of the in-plane normal stresses in the x and y direction, respectively, called and , and the in-plane shear stress . In Figure 5, the stress states of the bottom and top surfaces are visualized for a total load of 1000 N and 1800 N. These , and stress states are used to calculate the normal stress of each flaw for any load on the structure, as demonstrated in Equation (3). Combining this stress evolution with the corresponding flaw depths and the geometry factor as calculated in Equation (6), the evolution of the stress intensity factor for every flaw can be calculated and visualized, as depicted in Figure 6. The evolution of the stress intensity factor will show the location where a first flaw reaches a stress intensity factor equal to the fracture toughness of 0.75 MPa√m. Once this failing flaw is found, one knows the fracture load and the fracture location for this particular glass pane in a four-point bending setup. Figure 7 visualizes the fracture location for the random glass pane shown in the former figures. In this example, the fracture originated at the bottom surface location with and coordinates equal to 175 mm and 35 mm, respectively. These σ x , σ y and τ xy stress states are used to calculate the normal stress of each flaw for any load on the structure, as demonstrated in Equation (3). Combining this stress evolution with the corresponding flaw depths a i and the geometry factor Y as calculated in Equation (6), the evolution of the stress intensity factor K I for every flaw can be calculated and visualized, as depicted in Figure 6. The evolution of the stress intensity factor K I will show the location where a first flaw reaches a stress intensity factor equal to the fracture toughness of 0.75 MPa √ m. Once this failing flaw is found, one knows the fracture load and the fracture location for this particular glass pane in a four-point bending setup. Figure 7 visualizes the fracture location for the random glass pane shown in the former figures. In this example, the fracture originated at the bottom surface location with x and y coordinates equal to 175 mm and 35 mm, respectively.

Probability Density Distribution
Once the critical flaw is found, the fracture load and location are known. The crack calculation for a random glass pane can be repeated multiple times to obtain a histogram

Probability Density Distribution
Once the critical flaw is found, the fracture load and location are known. The crack calculation for a random glass pane can be repeated multiple times to obtain a histogram

Probability Density Distribution
Once the critical flaw is found, the fracture load and location are known. The crack calculation for a random glass pane can be repeated multiple times to obtain a histogram of the outcomes. Figure 8 presents the normalized histogram of this fracture load for 5000 repetitions, with a flaw density of two flaws for each square centimeter. The x-axis represents the fracture load, and the y-axis represents the observed probability density. The latter is a normalized value calculated by Equation (10), where w i is the bins' width in the distribution, c i represents the number of observations in the respective bin, and N is the total number of tests. The corresponding fracture locations for 5000 repetitions are illustrated in Figure 9.

Probability Density Distribution
Once the critical flaw is found, the fracture load and location are known. The crack calculation for a random glass pane can be repeated multiple times to obtain a histogram of the outcomes. Figure 8 presents the normalized histogram of this fracture load for 5000 repetitions, with a flaw density of two flaws for each square centimeter. The x-axis represents the fracture load, and the y-axis represents the observed probability density. The latter is a normalized value calculated by Equation (10), where is the bins' width in the distribution, represents the number of observations in the respective bin, and is the total number of tests. The corresponding fracture locations for 5000 repetitions are illustrated in Figure 9.
• (10)  To validate this outcome, the numerical results are compared to experimental results in Figure 8. Osnes et al. [13] performed 30 tests in the presented four-point bending setup and recorded the applied load at failure for each test. Comparing the numerical and

Probability Density Distribution
Once the critical flaw is found, the fracture load and location are known. The crack calculation for a random glass pane can be repeated multiple times to obtain a histogram of the outcomes. Figure 8 presents the normalized histogram of this fracture load for 5000 repetitions, with a flaw density of two flaws for each square centimeter. The x-axis represents the fracture load, and the y-axis represents the observed probability density. The latter is a normalized value calculated by Equation (10), where is the bins' width in the distribution, represents the number of observations in the respective bin, and is the total number of tests. The corresponding fracture locations for 5000 repetitions are illustrated in Figure 9.
• (10)  To validate this outcome, the numerical results are compared to experimental results in Figure 8. Osnes et al. [13] performed 30 tests in the presented four-point bending setup and recorded the applied load at failure for each test. Comparing the numerical and To validate this outcome, the numerical results are compared to experimental results in Figure 8. Osnes et al. [13] performed 30 tests in the presented four-point bending setup and recorded the applied load at failure for each test. Comparing the numerical and experimental results, it can be concluded that the numerical model described in this paper predicts the strength of glass panes properly. However, the experimental results indicate a larger frequency of small fracture loads and overall show a larger variability. At this point, there is no clear justification for this. Furthermore, Osnes et al. obtained a similar discrepancy with the experimental data. This is thus a point requiring further research. With this information, the description of the strength prediction model is considered completed.
Next, another parameter can be derived from the numerical model: the stress perpendicular to the flaw orientation of the failing flaw. This normal stress is an important parameter since it indicates the tensile capacity of the subjected flaws. The observed probability density and the cumulative density of this normal stress are shown in Figure 10. This normal stress distribution will be used to perform a convergence study and several parameter studies described later in this paper.
ancy with the experimental data. This is thus a point requiring further research. With this information, the description of the strength prediction model is considered completed.
Next, another parameter can be derived from the numerical model: the stress perpendicular to the flaw orientation of the failing flaw. This normal stress is an important parameter since it indicates the tensile capacity of the subjected flaws. The observed probability density and the cumulative density of this normal stress are shown in Figure 10. This normal stress distribution will be used to perform a convergence study and several parameter studies described later in this paper.

Convergence Study
To ensure that the strength prediction model's output is not affected by the number of simulations, a convergence study is performed. In this study, the output of the numerical model will converge to a repeatable solution when increasing the number of simulations. As mentioned, the normal stress distribution will be used to perform this convergence study. Table 2 lists the models used for this convergence study. For each model, the number of simulations is mentioned together with the calculated mean value of the normal stress distribution, the standard deviation and the coefficient of variation (COV). Figure 11 visualizes the convergence of the mean value and the standard deviation. From Table 2 and Figure 11, from 5000 simulations onwards, a reliable result in the overall shape of the distribution and its corresponding mean value and standard deviation can be obtained.

Convergence Study
To ensure that the strength prediction model's output is not affected by the number of simulations, a convergence study is performed. In this study, the output of the numerical model will converge to a repeatable solution when increasing the number of simulations. As mentioned, the normal stress distribution will be used to perform this convergence study. Table 2 lists the models used for this convergence study. For each model, the number of simulations is mentioned together with the calculated mean value of the normal stress distribution, the standard deviation and the coefficient of variation (COV). Figure 11 visualizes the convergence of the mean value and the standard deviation. From Table 2 and Figure 11, from 5000 simulations onwards, a reliable result in the overall shape of the distribution and its corresponding mean value and standard deviation can be obtained. (a) (b) Figure 11. Convergence of (a) the mean value and (b) the standard deviation. Figure 12 presents the PDFs and the CDFs of the above-mentioned models with the number of simulations starting from 5000 simulations. This figure clearly illustrates that even more simulations are needed to obtain a reliable result at the lower quantiles. To calculate the required number of simulations for a specific quantile, Equation (11) is adopted, and the corresponding results are displayed in Table 3. In Equation (11), N represents the number of simulations needed, is the quantile of interest and is the Figure 11. Convergence of (a) the mean value and (b) the standard deviation.  Figure 12 presents the PDFs and the CDFs of the above-mentioned models with the number of simulations starting from 5000 simulations. This figure clearly illustrates that even more simulations are needed to obtain a reliable result at the lower quantiles. To calculate the required number of simulations for a specific quantile, Equation (11) is adopted, and the corresponding results are displayed in Table 3. In Equation (11), N represents the number of simulations needed, P q is the quantile of interest and COV q is the target coefficient of variation for the quantile probability (set equal to 0.14 here to obtain a comparable number of simulations as in Figure 12). Such a quantile evaluation is relevant when the theoretical distribution describing the glass strength is unknown. However, as soon as a theoretical distribution is assigned, quantiles can be easily evaluated as soon as the parameters of the distribution (e.g., mean and standard deviation) are known.
(a) (b) Figure 11. Convergence of (a) the mean value and (b) the standard deviation. Figure 12 presents the PDFs and the CDFs of the above-mentioned models with the number of simulations starting from 5000 simulations. This figure clearly illustrates that even more simulations are needed to obtain a reliable result at the lower quantiles. To calculate the required number of simulations for a specific quantile, Equation (11) is adopted, and the corresponding results are displayed in Table 3. In Equation (11), N represents the number of simulations needed, is the quantile of interest and is the target coefficient of variation for the quantile probability (set equal to 0.14 here to obtain a comparable number of simulations as in Figure 12). Such a quantile evaluation is relevant when the theoretical distribution describing the glass strength is unknown. However, as soon as a theoretical distribution is assigned, quantiles can be easily evaluated as soon as the parameters of the distribution (e.g., mean and standard deviation) are known.

Model Selection and Distribution Fitting-AIC Calculation
To generalize the results, a distribution must be plotted on the outcome of the SPM, i.e., the probability density and the cumulative density of the normal stresses. Osnes [13] plotted a normal distribution of the data outcome of the SPM. It is observed that this normal distribution fits well with the data of a four-point bending setup with a length of 300 mm and a width of 60 mm. However, the data of a four-point bending setup with a length of 100 mm and a width of 20 mm differ from the normal distribution fit. This highlights the importance of investigating the right theoretical model and its possible dependence on the simulated case.
To select the right model, several criteria are available such as the Bayesian information criterion (BIC), the Akaike information criterion (AIC), and the minimum description length (MDL) [15]. In this paper, the Akaike information criterion (AIC) is chosen since model estimation and model selection are simultaneously accomplished [15]. The criterion was introduced in 1971 by Hirotugu Akaike, and the criterion is defined by an AIC value calculated according to Equation (12) for each model [25]. In Equation (12), k is the number of estimated parameters in the model andL is the maximum value of the likelihood function for the model. Once the AIC values of the competing models are calculated, the most appropriate model is that with the lowest AIC value.
In this paper, two-parameter distributions will be evaluated, which makes the value k equal to 2. The competing models are a normal distribution, a lognormal distribution, a gamma distribution, a Gumbel distribution and a two-parameter Weibull distribution.

Results and Discussion
The information in Section 2.3 allows us to select the most appropriate model to represent the data of the discussed case study in Section 2.2.3. Figure 13 indicates the normalized AIC values for the competing distributions, as discussed in Section 2.3, fitted on the data of Figure 10. From Figure 13, it can be concluded that the gamma distribution has the lowest AIC value, resulting in the most appropriate model to represent the data of Figure 10. Figure 14 illustrates the data of Figure 10 together with the five different models. It is clearly visible that the gamma and lognormal distributions are the best fit, especially for the lower quantiles. The normal distribution underestimates the strength (normal stress at fracture) for the lower quantiles.
the importance of investigating the right theoretical model and its possible dependence on the simulated case.
To select the right model, several criteria are available such as the Bayesian information criterion (BIC), the Akaike information criterion (AIC), and the minimum description length (MDL) [15]. In this paper, the Akaike information criterion (AIC) is chosen since model estimation and model selection are simultaneously accomplished [15]. The criterion was introduced in 1971 by Hirotugu Akaike, and the criterion is defined by an AIC value calculated according to Equation (12) for each model [25]. In Equation (12), k is the number of estimated parameters in the model and is the maximum value of the likelihood function for the model. Once the AIC values of the competing models are calculated, the most appropriate model is that with the lowest AIC value.
In this paper, two-parameter distributions will be evaluated, which makes the value equal to 2. The competing models are a normal distribution, a lognormal distribution, a gamma distribution, a Gumbel distribution and a two-parameter Weibull distribution.

Results and Discussion
The information in Section 2.3 allows us to select the most appropriate model to represent the data of the discussed case study in Section 2.2.3. Figure 13 indicates the normalized AIC values for the competing distributions, as discussed in Section 2.3, fitted on the data of Figure 10. From Figure 13, it can be concluded that the gamma distribution has the lowest AIC value, resulting in the most appropriate model to represent the data of Figure  10. Figure 14 illustrates the data of Figure 10 together with the five different models. It is clearly visible that the gamma and lognormal distributions are the best fit, especially for the lower quantiles. The normal distribution underestimates the strength (normal stress at fracture) for the lower quantiles.  With the validated strength prediction model and the above-mentioned information regarding the convergence and distribution fitting, different parameter studies are performed to give an idea of the most sensitive parameters in the model. The parameters that will be examined below are, in order of appearance, the scale factor, the ratio width/loading span, the ratio loading span/support span, the flaw density, the flaw shape, the maximum flaw depth and the fracture toughness. With the validated strength prediction model and the above-mentioned information regarding the convergence and distribution fitting, different parameter studies are performed to give an idea of the most sensitive parameters in the model. The parameters that will be examined below are, in order of appearance, the scale factor, the ratio width/loading span, the ratio loading span/support span, the flaw density, the flaw shape, the maximum flaw depth and the fracture toughness.
For every parameter, a table will present the varying parameter value together with the corresponding mean value, standard deviation, and COV for the normal stress at failure and the best fit (i.e., the distribution with the lowest AIC value). The setup introduced in Section 2 serves as a reference. This reference is displayed in every table of the parameter studies by means of a table footer.

Scale Factor
The investigation is a study on the total scaling of the setup. The width and the support span are scaled with the same factor, while the loading span is kept constant as half of the support span. The obtained results are presented in Table 4, the corresponding PDFs and CDFs are visualized in Figure 15, and the trends of the normalized AIC values, the mean values and the standard deviations are presented in Figure 16.  With the validated strength prediction model and the above-mentioned information regarding the convergence and distribution fitting, different parameter studies are performed to give an idea of the most sensitive parameters in the model. The parameters that will be examined below are, in order of appearance, the scale factor, the ratio width/loading span, the ratio loading span/support span, the flaw density, the flaw shape, the maximum flaw depth and the fracture toughness.
For every parameter, a table will present the varying parameter value together with the corresponding mean value, standard deviation, and COV for the normal stress at failure and the best fit (i.e., the distribution with the lowest AIC value). The setup introduced in Section 2 serves as a reference. This reference is displayed in every table of the parameter studies by means of a table footer.

Scale Factor
The investigation is a study on the total scaling of the setup. The width and the support span are scaled with the same factor, while the loading span is kept constant as half of the support span. The obtained results are presented in Table 4, the corresponding PDFs and CDFs are visualized in Figure 15, and the trends of the normalized AIC values, the mean values and the standard deviations are presented in Figure 16.    Figure 16b conclude that the mean value, the standard deviation, and COV decrease when a larger setup is considered, and vice versa. The larger the setup, the more flaws are present and are subjected to large tensile stresses in the loading span. This ensures a larger probability for a deep flaw to be present, which causes the glass to fail at a lower fracture load, which goes hand in hand with lower normal stresses at the moment of failure, as was expected. This trend is also clearly visible in Figure 15; the higher the scaling factor, the more the distribution shifts to the left and becomes more narrow.
Another observation is the shift in the best distribution fit. From Table 4, it can be concluded that the data of the small setup (i.e., the setup with a scaling factor of 0.25) fits best to a Gumbel distribution. When looking at the data of the setups with a scaling factor of 0.5-1, the best fit changes to a gamma distribution, and for the data of the biggest setups (i.e., the setups with a scaling factor of 2-4), to a normal distribution. This trend is also noticeable in Figure 16a, where all the normalized AIC values of the different models are summarized. This figure highlights that the Gumbel fit gets worse when approaching larger setups. The Weibull fit is the other way around; this fit gets better when approaching larger setups. For the performed scale factors, the gamma and lognormal fits perform well overall.    Figure 16b conclude that the mean value, the standard deviation, and COV decrease when a larger setup is considered, and vice versa. The larger the setup, the more flaws are present and are subjected to large tensile stresses in the loading span. This ensures a larger probability for a deep flaw to be present, which causes the glass to fail at a lower fracture load, which goes hand in hand with lower normal stresses at the moment of failure, as was expected. This trend is also clearly visible in Figure 15; the higher the scaling factor, the more the distribution shifts to the left and becomes more narrow.
Another observation is the shift in the best distribution fit. From Table 4, it can be concluded that the data of the small setup (i.e., the setup with a scaling factor of 0.25) fits best to a Gumbel distribution. When looking at the data of the setups with a scaling factor of 0.5-1, the best fit changes to a gamma distribution, and for the data of the biggest setups (i.e., the setups with a scaling factor of 2-4), to a normal distribution. This trend is also noticeable in Figure 16a, where all the normalized AIC values of the different models are summarized. This figure highlights that the Gumbel fit gets worse when approaching larger setups. The Weibull fit is the other way around; this fit gets better when approaching larger setups. For the performed scale factors, the gamma and lognormal fits perform well overall.

Width/Loading Span
The next parameter is the ratio of the width to the loading span. The width and the loading span in the reference setup are equal to 60 mm and 140 mm, respectively, resulting in a ratio of 0.43. The other models are made while keeping the loading span constant (i.e., 140 mm) and varying the width parameter to influence the ratio. Table 5

Width/Loading Span
The next parameter is the ratio of the width to the loading span. The width and the loading span in the reference setup are equal to 60 mm and 140 mm, respectively, resulting in a ratio of 0.43. The other models are made while keeping the loading span constant (i.e., 140 mm) and varying the width parameter to influence the ratio. Table 5 displays the data of the different models, the corresponding PDFs and CDFs are visualized together in Figure 17, and the trends of the normalized AIC values, the mean values and the standard deviations are presented in Figure 18.  The table and Figure 18b display a decrease in the mean value, standard deviation and COV when increasing the ratio, and vice versa. For these simulations, a larger ratio is caused by a larger width, resulting in a larger region subjected to the maximum bending moment. This means that more flaws are present in the zone with the maximum stresses, increasing the probability of the glass failing at lower stresses. This trend is also clearly visible in Figure 17; the higher the ratio, the more the distribution shifts to the left and becomes more narrow.
Observing the column with the best fit in Table 5, a shift for the best fit is noticeable. For the small ratios (i.e., slender samples), a lognormal distribution has the lowest AIC value. Ratio 0.2 to 0.4 corresponds best to a gamma distribution. In the range from 0.5 to 0.7, the AIC values for a normal fit and a gamma fit are very close to each other, so both fits can be used. The largest ratios, starting from 0.8, have the normal distribution as the lowest AIC value. Figure 18a illustrates the same trend. Additionally, the normalized AIC values visualize that the Gumbel fit becomes worse when approaching larger values for the ratio of the width to the loading span. Again, the gamma and lognormal fit perform well for all the parameter values considered. The table and Figure 18b display a decrease in the mean value, standard deviation and COV when increasing the ratio, and vice versa. For these simulations, a larger ratio is caused by a larger width, resulting in a larger region subjected to the maximum bending moment. This means that more flaws are present in the zone with the maximum stresses, increasing the probability of the glass failing at lower stresses. This trend is also clearly visible in Figure 17; the higher the ratio, the more the distribution shifts to the left and becomes more narrow.
Observing the column with the best fit in Table 5, a shift for the best fit is noticeable. For the small ratios (i.e., slender samples), a lognormal distribution has the lowest AIC value. Ratio 0.2 to 0.4 corresponds best to a gamma distribution. In the range from 0.5 to 0.7, the AIC values for a normal fit and a gamma fit are very close to each other, so both fits can be used. The largest ratios, starting from 0.8, have the normal distribution as the lowest AIC value. Figure 18a illustrates the same trend. Additionally, the normalized AIC values visualize that the Gumbel fit becomes worse when approaching larger values for the ratio of the width to the loading span. Again, the gamma and lognormal fit perform well for all the parameter values considered. The table and Figure 18b display a decrease in the mean value, standard deviation and COV when increasing the ratio, and vice versa. For these simulations, a larger ratio is caused by a larger width, resulting in a larger region subjected to the maximum bending moment. This means that more flaws are present in the zone with the maximum stresses, increasing the probability of the glass failing at lower stresses. This trend is also clearly visible in Figure 17; the higher the ratio, the more the distribution shifts to the left and becomes more narrow.
Observing the column with the best fit in Table 5, a shift for the best fit is noticeable. For the small ratios (i.e., slender samples), a lognormal distribution has the lowest AIC value. Ratio 0.2 to 0.4 corresponds best to a gamma distribution. In the range from 0.5 to 0.7, the AIC values for a normal fit and a gamma fit are very close to each other, so both fits can be used. The largest ratios, starting from 0.8, have the normal distribution as the lowest AIC value. Figure 18a illustrates the same trend. Additionally, the normalized AIC values visualize that the Gumbel fit becomes worse when approaching larger values for the ratio of the width to the loading span. Again, the gamma and lognormal fit perform well for all the parameter values considered.

Loading Span/Support Span
Another parameter study is performed on the ratio of the loading span to the support span. The loading and support span in the reference setup is equal to 140 mm and 280 mm, resulting in a ratio of 0.5. The other models are made while keeping the support span constant (i.e., 280 mm) and varying the loading span parameter to influence the ratio. Table 6 displays the data for the models, all PDFs and all CDFs are visualized together in Figure 19, and the trends of the normalized AIC values, the mean values and the standard deviations are presented in Figure 20.

Loading Span/Support Span
Another parameter study is performed on the ratio of the loading span to the support span. The loading and support span in the reference setup is equal to 140 mm and 280 mm, resulting in a ratio of 0.5. The other models are made while keeping the support span constant (i.e., 280 mm) and varying the loading span parameter to influence the ratio. Table 6 displays the data for the models, all PDFs and all CDFs are visualized together in Figure 19, and the trends of the normalized AIC values, the mean values and the standard deviations are presented in Figure 20.    Figure 20b indicate a smaller mean value, standard deviation and COV for larger ratios, and vice versa. The larger the loading span, the larger the region subjected to the maximum stresses. This way, the number of flaws in this region is higher, increasing the probability of the glass failing at smaller stresses. Again, this trend is also clearly visible in Figure 19; the larger the ratio, the more the distribution shifts to the left and becomes more narrow.   Figure 20b indicate a smaller mean value, standard deviation and COV for larger ratios, and vice versa. The larger the loading span, the larger the region subjected to the maximum stresses. This way, the number of flaws in this region is higher, increasing the probability of the glass failing at smaller stresses. Again, this trend is also clearly visible in Figure 19; the larger the ratio, the more the distribution shifts to the left and becomes more narrow.
A study of the column with the best fit tells us a shift in the best distribution fit. The data of the ratios below and equal to 0.5 all list the best fit with a gamma distribution. From the data of the ratios higher and equal to 0.6, both a normal distribution and a gamma distribution are good fits since the AIC values hardly differ. Figure 20a illustrates the same trend, with the extra information that the Gumbel distribution worsens when considering higher ratios. Overall, the gamma fit performs well.

Flaw Density
In Section 2, the flaw density was presented as an integer. However, from the literature, it is clear that this parameter varies considerably, depending on the cutting and finishing processes used. Wereszczak et al., (2014) [22] investigated the flaw density for soda-limesilica glass when scored and bent and cut with a water jet. The authors found that the flaw density can vary from 1.18 to 2.60 flaws/cm 2 . [22] This shows the importance of a parameter study on the flaw density. Table 7 displays the models used for this parameter study, all PDFs and all CDFs are visualized together in Figure 21, and Figure 22 presents the trends of the normalized AIC values, the mean values and the standard deviations.      Figure 22b prove that a larger flaw density causes a smaller mean value, standard deviation and COV for the normal stress distribution. A larger flaw density resulted in more flaws on the glass surface and subjected to the tensile stresses resulting from the bending moments. This ensures that there is a larger probability for a deep flaw to be present, which causes the glass to fail at a lower fracture load, which goes hand in hand with lower normal stresses at the moment of failure. Figure 21 visualizes the smaller values and a more narrow distribution for larger flaw densities.
Another observation is again a switch in the best distribution fit. The data for a setup with a flaw density equal to 0.5 flaws/cm 2 has the best fit with a lognormal distribution. For flaw densities ranging from 1.0 to 2.5 flaw/cm 2 , the best fit changes to a gamma distribution. Furthermore, flaw densities higher than and equal to 3.0 flaws/cm 2 have a best fit with a normal distribution. The AIC values in Figure 22a also illustrate the Gumbel fit worsening when considering larger flaw densities.   Figure 22b prove that a larger flaw density causes a smaller mean value, standard deviation and COV for the normal stress distribution. A larger flaw density resulted in more flaws on the glass surface and subjected to the tensile stresses resulting from the bending moments. This ensures that there is a larger probability for a deep flaw to be present, which causes the glass to fail at a lower fracture load, which goes hand in hand with lower normal stresses at the moment of failure. Figure 21 visualizes the smaller values and a more narrow distribution for larger flaw densities.
Another observation is again a switch in the best distribution fit. The data for a setup with a flaw density equal to 0.5 flaws/cm 2 has the best fit with a lognormal distribution. For flaw densities ranging from 1.0 to 2.5 flaw/cm 2 , the best fit changes to a gamma distribution. Furthermore, flaw densities higher than and equal to 3.0 flaws/cm 2 have a best fit with a normal distribution. The AIC values in Figure 22a also illustrate the Gumbel fit worsening when considering larger flaw densities.

Flaw Shape
In the Section above, the flaw shape a/c is taken equal to 1 based on the work of Levengood [23]. Consulting Equations (6)-(9), as presented in Section 2.2.2, this assumption leads to a geometry factor of 0.729. However, this is an assumption where all flaws are idealized as elliptic. In reality, much more flaws will be irregularly shaped [12]. Therefore, in this parameter study, the flaw shape is varied to investigate the influence of more narrow and elongated flaws. Table 8 displays the models used for this parameter study; the respective PDFs and CDFs are illustrated in Figure 23. Figure 24 presents the trends of the normalized AIC values, the mean values and the standard deviations. tion leads to a geometry factor of 0.729. However, this is an assumption where all flaws are idealized as elliptic. In reality, much more flaws will be irregularly shaped [12]. Therefore, in this parameter study, the flaw shape is varied to investigate the influence of more narrow and elongated flaws. Table 8 displays the models used for this parameter study; the respective PDFs and CDFs are illustrated in Figure 23. Figure 24 presents the trends of the normalized AIC values, the mean values and the standard deviations.  Levengood [23]. Consulting Equations (6)-(9), as presented in Section 2.2.2, this assumption leads to a geometry factor of 0.729. However, this is an assumption where all flaws are idealized as elliptic. In reality, much more flaws will be irregularly shaped [12]. Therefore, in this parameter study, the flaw shape is varied to investigate the influence of more narrow and elongated flaws. Table 8 displays the models used for this parameter study; the respective PDFs and CDFs are illustrated in Figure 23. Figure Table 8 clearly demonstrates a minimum in the mean value, the standard deviation and the COV of the normal stress distribution when the flaw shape equals 0.8. This trend is extra visible in Figure 24b, where a plot is made of the mean values of normal stress at failure in function of the flaw shape. It results from Equations (6)- (9). These Equations calculate the corresponding geometry factor for the given flaw shape, and a maximum geometry factor of 0.734 is reached for a flaw shape equal to 0.833. Looking at Equation (2), this maximum leads to a higher stress intensity factor which means it will reach the fracture toughness sooner, resulting in a lower fracture load and, thus lower normal stresses at failure (i.e., a minimum for a flaw shape of 0.833).
The AIC calculation illustrates that there is no shift in the best fit for this parameter study. A normal distribution and a gamma distribution are good fits for all the models in this parameter study.

Maximum Flaw Depth
The maximum flaw depth a max is the parameter used to simulate the depths of all the flaws that must be distributed over the glass surfaces. This parameter equals 0.1 mm, but the literature has already demonstrated a large spread for this value in reality. In the paper of Wereszczak et al. [22], different scanned glass surfaces showed the largest identified flaw on a glass pane going from 105 µm to 398 µm. This demonstrates the importance of performing a parameter study on this maximum flaw depth parameter. Table 9 displays the models used for this parameter study; all PDFs and all CDFs are visualized together in Figure 25. The trends of the normalized AIC values, the mean values and the standard deviations are presented in Figure 26. It can be concluded that the mean value and the standard deviation of the normal stress distribution decrease when the maximum flaw depth is increasing, and vice versa. The obtained COV is, however, approximately constant. Referring to Equation (2), a larger flaw depth results in a larger stress intensity factor which means it will reach the fracture toughness sooner, resulting in a lower fracture load and, thus lower normal stresses.  The AIC calculation again shows no shift in the best fit for this parameter study. Both a normal distribution and a gamma distribution are again good fits for all the models in this parameter study. Table 9. Data for a varying maximum flaw depth. The AIC calculation again shows no shift in the best fit for this parameter study. Both a normal distribution and a gamma distribution are again good fits for all the models in this parameter study. Table 9. Data for a varying maximum flaw depth. The AIC calculation again shows no shift in the best fit for this parameter study. Both a normal distribution and a gamma distribution are again good fits for all the models in this parameter study.

Fracture Toughness
The last parameter study is the one on fracture toughness. The literature has proven that, at inert conditions, a fixed value of 0.75 MPa √ m for the fracture toughness is a good assumption for soda-lime-silica glass [16,26,27]. Yet, different investigations have demonstrated that there is a certain variation in this parameter. Gong et al., (2001) [28] performed a statistical analysis of the fracture toughness of soda-lime-silica glass, determined by indentation. The experimental research displays a spread of the fracture toughness from 0.5 to 0.9 MPa √ m. [28]. The models used in this parameter study are based on the values of Gong et al., (2001). The results are displayed in Table 10, Figures 27 and 28.

Fracture Toughness
The last parameter study is the one on fracture toughness. The literature has proven that, at inert conditions, a fixed value of 0.75 MPa√m for the fracture toughness is a good assumption for soda-lime-silica glass [16,26,27]. Yet, different investigations have demonstrated that there is a certain variation in this parameter. Gong et al. (2001) [28] performed a statistical analysis of the fracture toughness of soda-lime-silica glass, determined by indentation. The experimental research displays a spread of the fracture toughness from 0.5 to 0.9 MPa√m. [28]. The models used in this parameter study are based on the values of Gong et al. (2001). The results are displayed in Table 10 The results demonstrate an increase in the mean value and the standard deviation of the data of the normal stresses when increasing the fracture toughness. Note that the COV is more or less constant, resulting in a linear increase for both the mean value and the standard deviation. The increase in the fracture toughness allows for larger stress intensity factors before the glass breaks. This results in larger allowable normal stresses.
The best fits do not experience a shift for varying fracture toughness. The fracture toughness is a deterministic value and does not influence the best fit. Any change is the

Fracture Toughness
The last parameter study is the one on fracture toughness. The literature has proven that, at inert conditions, a fixed value of 0.75 MPa√m for the fracture toughness is a good assumption for soda-lime-silica glass [16,26,27]. Yet, different investigations have demonstrated that there is a certain variation in this parameter. Gong et al. (2001) [28] performed a statistical analysis of the fracture toughness of soda-lime-silica glass, determined by indentation. The experimental research displays a spread of the fracture toughness from 0.5 to 0.9 MPa√m. [28]. The models used in this parameter study are based on the values of Gong et al. (2001). The results are displayed in Table 10 The results demonstrate an increase in the mean value and the standard deviation of the data of the normal stresses when increasing the fracture toughness. Note that the COV is more or less constant, resulting in a linear increase for both the mean value and the standard deviation. The increase in the fracture toughness allows for larger stress intensity factors before the glass breaks. This results in larger allowable normal stresses.
The best fits do not experience a shift for varying fracture toughness. The fracture toughness is a deterministic value and does not influence the best fit. Any change is the The results demonstrate an increase in the mean value and the standard deviation of the data of the normal stresses when increasing the fracture toughness. Note that the COV is more or less constant, resulting in a linear increase for both the mean value and the standard deviation. The increase in the fracture toughness allows for larger stress intensity factors before the glass breaks. This results in larger allowable normal stresses.
The best fits do not experience a shift for varying fracture toughness. The fracture toughness is a deterministic value and does not influence the best fit. Any change is the result of randomness in the Monte Carlo simulations. Again, both a normal distribution and a gamma distribution are good fits for all the models in this parameter study.

Conclusions
The strength prediction model by Osnes has been expanded with an AIC calculation to perform a model selection for every outcome of the SPM. This model selection proved that the gamma distribution was the best distribution fit for the data of the presented case study.
Based on a parameter study, it is concluded that the most appropriate model to represent the data of the normal stresses at failure was largely influenced by the following four parameters: the scale factor, the ratio of width to loading span, the ratio of loading span to support span, and the flaw density. As a general conclusion of this parametric study, the more flaws are subjected to the maximum tensile stresses, the more the Weibull distribution and the normal distribution become appropriate. Conversely, the less flaws are subjected to the maximum tensile stresses, the more the Gumbel distribution becomes the best fit. Overall, the gamma distribution performs well for all the simulated cases.
The parameter study further indicates that the mean value, standard deviation and COV of the glass panel strength increase while decreasing the scale factor, the ratio of the width to the loading span, the ratio of the loading span to the support span, or the flaw density. The mean value and the standard deviation of the data of the normal stresses increase more or less linearly, with a decrease in maximum flaw depth or an increase in fracture toughness. In these cases, the COV stays more or less constant. As a last conclusion, it was found that the mean value, the standard deviation and the COV of the data of the normal stresses reach a minimum when the scale factor of the flaw approaches a value of 0.833.
The performed parameter study and model selection gives a view of how the glass strength distribution changes while considering other parameter values. For example, the strength of a glass panel fresh from production (i.e., with a smaller flaw density) is more describable with a Gumbel distribution. In contrast, a glass panel already in service for several years (i.e., with a larger flaw density) corresponds more closely to a Weibull or normal distribution.
On the other hand, the performed parameter study of the flaw shape, the maximum flaw depth, and the fracture toughness highlight the importance of selecting the right values for the considered parameters. Until now, values have been chosen based on experimental results from the literature. However, this study emphasizes that a small difference in value can reasonably affect the strength distribution.
In further research, the flaw shape, the maximum flaw depth, and the fracture toughness should be implemented as probabilistic values to account for the uncertainty of these parameters. Additionally, other cases (e.g., coaxial double-ring setups) could be checked to investigate whether the same trends can be found.

Conflicts of Interest:
The authors declare no conflict of interest.