A numerical method for analysis of fracture statistics of glass and simulations of a double ring bending test

The results from a new numerical method for simulating the strength and fracture locations of small glass specimens subjected to double ring bending are compared with experimental data. The method implements the weakest-link principle while assuming the existence of Grifﬁth ﬂaws. A Weibull distribution for the strength is simulated based on a single popula-tionofParetodistributedcracksizes.Theeffectofusing different fracture criteria is investigated. An alternative distribution is simulated based on two populations of ﬂaws. This distribution models the apparent bimodality in the empirical data set. The numerical method is dependent on a representation of the surface ﬂaws condition in glass. As new techniques become available for examining the surface characteristics, this numerical method is promising as a means for modelling the strength better than current methods do.


Introduction
In order to explain and predict the strength of annealed glass a range of concepts and methods have been applied with mixed results. Typically, the strength is explained assuming the existence of Griffith flaws and supposing the weakest-link principle. Predictions are based either on some standard distribution or on tables and diagrams obtained using a modelling tool such as the Glass Failure Prediction Model (GFPM) (Beason and Morgan 1984). There is disagreement among researchers as to which prediction model is the correct one to use (Fischer-Cripps and Collins 1995). A range of experiments have shown a consistent bilinearity in the probability plots when the Weibull distribution is used for modelling the strength of annealed glass (Veer 2007;Veer et al. 2009). As regards the GFPM, it has been said that it "is best suited to representing glass strength for specific test conditions." (Reid 2007) Neither the standard distributions nor the GFPM are able to consistently provide for an acceptable goodness-of-fit while modelling data from experiments, something that is called for in a prediction model with true potential. At the same time, structural glass is gaining in popularity among designers and units are being installed in buildings and public spaces worldwide at an increasing rate. The search for a failure prediction model is therefore as topical as ever. Moreover, a study has indicated that shear stress might affect the observed strength of glass in double ring bending tests (Reid 2007). Shear stress is generally not considered in current failure models for glass. In this article, a numerical method for predicting the failure of glass is investigated and applied to double ring bending tests. The method is general and can be applied to a range of specimen geometries, loading setups and support conditions. The results are compared with experimental data. The presented method depends on a model that is based on fracture mechanics and the weakest-link principle while assuming a preexisting population of surface cracks. Stress corrosion is not directly considered in this study. The presented method should not be taken as a complete and final strength design tool. The aim is to promote a failure model for glass that is based on a logical and tractable representation of the surface condition in glass with a suitable consideration of the fracture mechanics. With further research, this could in the end lead to an improved strength design tool for use with glass.

Background
Flaws in glass are capable of promoting brittle failure due to the lack in capacity for plastic flow. While assuming that the surface contains a large number of minute flaws that act like cracks, so-called Griffith flaws, it is possible to explain the scatter in fracture location observed in experiments. It also helps to explain the variation in failure stress observed and the relatively low strength attained in practice. Surface flaws arise in the production line during manufacture as well as in subsequent handling, transportation, assembly, use, and maintenance. Bulk flaws are disregarded in the following, cf. Bourhis (2008). Griffith (1920) modelled crack growth as a reversible thermodynamical process. For a crack subjected to mode I opening displacement, fracture is governed by the following criterion where K I is the Stress Intensity Factor (SIF) and K Ic denotes the fracture toughness (Irwin 1957). The value of K Ic for glass has been estimated at about 0.75 MPa m 1/2 (Mencik 1992). It is assumed that the individual cracks do not interact with each other. The shape of a surface crack in glass is typically conceived of as being either a long, straight-fronted plane crack or a semi-circular crack (Haldimann 2006). There exist several solutions to the calculation of the SIF for a semicircular crack subjected to a uniform tensile stress field σ n oriented perpendicular to the crack plane. Accord-σ 1 σ 2 θ Fig. 1 A plane crack subjected to a biaxial stress field with the crack plane inclined at an angle θ in the coordinate system of the principal stresses σ 1 and σ 2 ing to one solution the SIF at the deepest point on the crack contour is (Newman and Raju 1981) where a denotes the crack depth, see also Thiemeier et al. (1991). Figure 1 illustrates a crack subjected to a biaxial stress field with the crack plane inclined at an angle θ in the coordinate system of the principal stresses σ 1 and σ 2 . If the crack plane is oriented perpendicular to the Maximum Principal Tensile Stress (MPTS) σ 1 then is substituted into Eq. (2). Otherwise, the tensile stress acting normal to the crack plane can be calculated as The presence of shear stress does not have any effect in a pure mode I fracture criterion. There exists a range of fracture criteria for a crack subjected to both normal and shear stresses while assuming mode I crack opening and mode II in-plane shearing displacements. One such mixed mode fracture criterion which is based on the maximum non-coplanar energy release rate (Hellen and Blackburn 1975) is given by the following inequality where the left-hand side of inequality (5) is a mode Iequivalent SIF, cf. Thiemeier et al. (1991). For a semicircular crack K II can be approximated as (Thiemeier et al. 1991) In Eq. (6) ν is Poisson's ratio. τ is the shear stress in the crack plane which can be calculated from the in-plane principal stresses as The idea that you can calculate the distribution of macroscopic strength of a stressed solid by starting from an analysis of the microscopic defects dates back at least to Peirce (1926). Peirce (1926) formulated the Weakest-Link Principle (WLP), i.e. that the strength of a chain is governed by its weakest link, and applied it in the study of the tensile strength of cotton yarn. Also using the WLP, Weibull (1939) came up with the following distribution function for the strength of a brittle solid where B, denoted "the risk of rupture", is a function of body size and tensile stress. According to Weibull (1939), a simple mathematical form that is in general accord with experimental data is where k and m denote the scale and shape parameters, respectively. Inserting Eq. (9) into Eq. (8) gives the standard two-parameter Weibull distribution function where k is also the 63rd percentile (Wachtman et al. 2009). Various derivations of the strength distribution function for a brittle solid are offered by e.g. Freudenthal (1968), Matthews et al. (1976), Batdorf and Heinisch (1978), Evans and Jones (1978) and Danzer (1992). In general, the derivation is based on a subdivision of the stressed solid into regions. It is assumed that there exists a population of nonoverlapping cracks which are distributed among the regions. Each crack is associated with a critical stress. It is assumed that the stress state varies slowly so that all cracks within a subdivided region are subjected to the same nominal stress. The solution methods, which are analytical, vary. Also varying are certain assumptions, such as whether or not it is supposed that the fracture of the crack depends only on the component of stress normal to the crack plane, whether or not there exist multiple crack populations, etc. Essentially, the analytical expression for the strength distribution is obtained through a limit operation in which the region size shrinks infinitesimally while the number of subdivided regions increases indefinitely. All these solution methods are capable of producing the fundamental Eq. (8). However, the mathematics soon become intractable when all but the simplest assumptions are made for the stress state, fracture criterion, crack size distribution, flaw density, crack plane orientation, and the existence of multiple flaw populations. Yankelevsky (2014) offers a numerical solution method to the problem of determining the strength distribution of a brittle solid while building upon the same general ideas as in the aforementioned studies except that the limit operation is not carried out. In other words, it is not necessary to assume that a crack of some finite size is contained within an infinitesimally small space, cf. Afferrante et al. (2006). Yankelevsky illustrates the method in a study of a glass square plate subjected to bending. He neglects bulk flaws and considers failures starting from the surface area only. The surface area of the plate is subdivided into unit cells measuring 1 cm 2 . One crack is distributed into each cell. The flaw size density function proposed by Yankelevsky (2014) and which is motivated for use with glass material can be interpreted as a truncated exponential distribution. The square plate is laterally supported along two opposite edges and subjected to a line-load at midspan producing a uniaxial state of stress in the plane of the tensioned surface. A Monte Carlo simulation is carried out for a large sample of thousands of virtual specimens. This numerical method offers a tractable way of calculating the strength distribution as well as the fracture location distribution for arbitrary stress states, fracture criteria, crack plane orientations, crack size distributions, and multiple flaw populations. However, in Yankelevsky (2014), only a uniaxial tensile stress field is considered where the cracks are stressed normal to their crack planes. Subcritical crack growth is not considered. Nor is the method applied to a double ring bending test which is quite a common and relatively inexpensive method to evaluate the strength of small glass plates (Dalgliesh and Taylor 1990). Based on Hertzian indentation tests it has been suggested that flaw size in glass can be closely fitted by a Pareto distribution (Poloniecki and Wilshaw 1971;Poloniecki 1974;Tandon et al. 2013). The Pareto distribution has the scale and shape parameters a 0 > 0 and c > 0 and the distribution function is (Forbes et al. 2010) It has moreover been shown that the Weibull distribution function is derived from the WLP if it is assumed that the surface flaws condition is represented by a single population of cracks whose size is Pareto distributed in the tail (Jayatilaka and Trustrum 1977). In this view the Weibull shape parameter is a true material parameter. Then, the relation between the shape parameters m and c of the Weibull and Pareto distributions, respectively, is found to be When stressed in an ambient atmosphere, glass strength is reduced over time due to a process known as static fatigue which is due to subcritical crack growth, the effects of which are only observed when the mode I SIF lies above a threshold limit value at about 0.25 MPa m 1/2 (Wiederhorn and Bolz 1970). In Charles' stress corrosion rate theory (Charles 1958a, b), subcritical crack growth is explained as a thermally activated chemical process whereby water moisture interacts with tensile stress at the crack tip. Equation (12), however, often approximates observed values of subcritical crack growth (Mencik 1992 where v is the subcritical crack growth velocity, A is a constant, and n is the stress corrosion parameter. While the value of n was repeatedly estimated at about 16 for soda-lime glass in ambient conditions, the value of A at 50% relative humidity was estimated in a range spanning more than two orders of magnitude, see Schula (2015) for an overview of those experiments. Hence, it is generally challenging to predict subcritical crack growth in ambient conditions.

Surface flaws concept
For the representation of the surface flaws condition, we consider two models. The first one comprises a single population of semi-circular edge cracks with a Pareto distributed crack size. The second model comprises a dual population of semi-circular edge cracks with a Pareto and Fréchet distributed crack size, respectively.
In both cases, a choice of crack density at 2 cm −2 is made. The purpose with the dual population model presented here is to provide a logical basis for a strength distribution with a bimodality. The choice of crack density at 2 cm −2 is guided by the following observation. Based on optical scanning techniques applied to a pair of small soda-lime silicate glass plates in the asreceived condition there were 632 flaws observed and it was noted that the flaw mean density varied between 1.2 and 2.6 cm −2 for flaw sizes greater than approximately 8 microns (Wereszczak et al. 2014).

Single population model
For the single population model, it is assumed that the cracks are uniformly distributed over the surface area of the original plate and that the crack planes are oriented between [0, π) according to a uniform distribution.
The logical basis for the selected choice of single population model are the Hertzian indentation tests that have been carried out in the past (Poloniecki and Wilshaw 1971;Poloniecki 1974;Tandon et al. 2013) and which have provided data that could be closely fitted by a Pareto distribution, see Sect. 2.

Dual population model
For the dual population model it is assumed that it comprises two populations of semi-circular edge cracks with a Pareto and Fréchet distributed crack size, respectively. All cracks are uniformly distributed over the surface area of the original plate and the crack planes are oriented perpendicular to the MPTS. The Pareto population cracks represent large surface flaws. The Fréchet cracks represent small surface flaws. It is assumed that the number of Pareto cracks is a small fraction of the total number of cracks. It is assumed that the fraction is 0.002.
The logical basis for the dual population model is the following. First, glass fracture statistics tend to produce bimodalities in the probability plots according to e.g. Veer et al. (2009). In fact, the experiment considered in Sect. 5 is no exception because the histogram of the data appears to exhibit two modes, see also Fig. 3. Other researchers have suggested to represent the surface cracks using two populations. Mencik (1992) distinguishes between several populations of surface flaws according to their origin. In doing so, Mencik (1992) distinguishes between a large flaws and a small flaws population of cracks with relevance for the practical engineering strength of glass. Mencik (1992) characterizes the large flaws population as being responsible for the tensile stress to decrease to 20-60 MPa. Substituting these values into Eq. (2) while assuming that the SIF equals to 0.75 MPa m 1/2 , the corresponding crack depth is found to be 94-850 microns. Mencik (1992) characterizes the small flaws population as containing cracks smaller than a hundredth of a millimeter in depth. He associates this with a strength reduced to 60-200 MPa. Substituting these values into Eq. (2) yields a corresponding crack depth of 8-94 microns. A statistical model for characterizing glass strength when two flaw populations are superimposed due to abrasive phenomena has been proposed in Pisano and Carfagni (2017). Pathirana et al. (2017) implemented a dual population of cracks in a numerical model for the evaluation of the strength distribution in panels subjected to point contact actions. Second, the choice for the value of the fraction of large cracks, i.e. 0.002, is guided by the following observation. Out of the total number of flaws detected in the investigations by Wereszczak et al. (2014), the proportion of large flaws greater in size than or equal to about 200 microns was approximately 0.002. This corresponds with a crack depth of 100 microns assuming that the flaws are semicircular surface cracks. Taking a crack depth of about 100 microns as a value that separates large flaws from small flaws is through adoption of the line of reasoning in Mencik (1992). Third, the logical basis for the Pareto distribution are the Hertzian indentation tests that have been carried out in the past and which have been mentioned earlier in this paper already. Fourth, the logical basis for the Fréchet distribution is motivated as follows. Assuming that small flaws are exceedingly numerous, one might select only the greatest small flaw in a given region and let this one determine the fracture mechanical behaviour of the small flaws population in that region (Freudenthal 1968). Because it is assumed that the cracks in the small population are abundant, it is supposed that the selected crack plane is oriented approximately normal to the maximum principal tensile stress. If the numerous small flaws have an independent Pareto size distribution then in the limit the largest flaw size is Fréchet distributed (Beirlant et al. 2004). For extreme-value theoretical reasons the greatest flaw size among a large set of flaws whose size is iden-tically and independently distributed is approximately Fréchet distributed if the following holds (Horst 2009); the sampled distribution has a range which is unlimited from above and its distribution function F is such that there exist some positive numbers k, A such that Hence, the Pareto distribution lies in the domain of attraction of the Fréchet distribution (Beirlant et al. 2004).

Numerical modelling tool
Here follows a description of a numerical modelling tool for the strength of glass plates in bending based on an implementation of the weakest-link principle and some concept for the surface flaws condition. The numerical method adopted in this study is based on the Monte Carlo simulation method carried out by Yankelevsky (2014). The most important difference between the present study and Yankelevsky (2014) is that the present study considers multiple flaw populations with arbitrary crack plane orientations and a mixed mode fracture criterion. Float glass is usually produced and shipped in a standard size so-called jumbo plate with the dimensions 3.21 × 6.00 m 2 . Taking the standard jumbo plate as a starting point, the plate is subdivided into unit cells of 1 mm 2 . This cell size provides a reasonable compromise between resolution and computational cost. A set of flaws are randomly scattered across the cells according to a uniform distribution, although in general another spatial distribution could be adopted. It is supposed that the stochastic orientation of the crack planes is uniformly distributed. This assumption might not be conservative, however, if there is a tendency for the flaws to lie in some particular direction due to e.g. machining abrasion or contact with the rollers during manufacture. The total number of flaws on the jumbo plate is fixed and depends on the flaw density. It is assumed that the flaw density is 2 cm −2 yielding a total of 385,200 flaws on either face of the plate. Each flaw is independently assigned a size based on some statistical distribution function which depends on the particular flaws concept that is adopted. The random flaws are resampled in each new simulation of the jumbo plate. The stochastics of the flaws are the location, the orientation in the plane, and the size. Next, a specimen of given dimensions is extracted and separated from the jumbo plate. In the following the cut out specimen is analysed.
The cut out specimen is subjected to an arbitrary loading in increments and the stress field history at the centre of each flaw-containing unit cell determined. It is supposed that the load type is such that tensile stress actually develops on the face of the cut out specimen otherwise failure will not be detected based on fracture mechanics. In general the stress response is non-linear. The loading increment is chosen so as to produce a tensile stress increase of no more than 1 MPa per increment anywhere on the specimen. However, if the response is linear then it suffices with one increment and to scale the results. The complete stress history needs only to be calculated once for a given loading type and specimen geometry because the stochastics of the cracks do not affect the distribution of nominal bending stresses. It is assumed that the sum of load increments is sufficiently large in relation to the given flaw characteristics, i.e. flaw density, flaw size distribution, etc., to prompt fracture. Otherwise, failure might not have been detected by the end of the last load increment. There exists a SIF envelope that meets with the fracture toughness at some point in time, the smallest of which is identified as the time of failure. If the crack planes are always oriented perpendicular to the MPTS then the SIF envelope is calculated using Eqs. (2) and (3). For reference, this case is denoted MPTS mode I fracture criterion. If the crack planes are inclined at an oblique angle in the coordinate system of the principal stresses while mode I opening displacement is considered then the SIF envelope is calculated using Eqs. (2) and (4). This case is denoted oblique angle mode I fracture criterion. If both mode I opening and mode II shearing displacements are accounted for then the SIF envelope is calculated using the left-hand side of inequality (5). This is the mixed mode fracture criterion. By token of the WLP, the fracture origin is determined from the first unit cell that contains a flaw with a SIF exceeding the fracture toughness. A search algorithm is used to detect this cell. By carrying out simulations on a whole series of cut out specimens it is possible to obtain a sample of the fracture stress which is defined as the MPTS at the failure origin. In this study the number of cut out specimens in a simulation series is 5,000. This sample size offers a reasonable compromise between precision and computational cost. Figure 2 illustrates the main steps leading up to the creation of the SIF envelope; (a) the flaw stochastics are sampled, (b) the in-plane principal stresses are determined at each load increment and (c) the SIF envelope is calculated per load increment based on either of the fracture criteria. Failure is prompted at the first instance of intersection between SIF envelope and fracture toughness (white disc). Likewise, the failure origin is determined by the first unit cell that contains a SIF which exceeds the fracture toughness. The so-called critical stress is the uniform tensile stress perpendicular to a given crack plane that would bring about failure with a pure mode I fracture criterion. The critical stress can be calculated with Eqs. (1) and (2).

Experimental data comparison
Double ring bending tests are frequently carried out to evaluate the strength of glass. In this testing device a glass plate is supported on a reaction ring and subjected to an applied loading through a smaller concentric ring on its opposite side. An equibiaxial state of stress is produced within the loading ring. Failures that start from edges are eliminated because tensile stress diminishes near the edges. Some experimenters discard any observation associated with a failure originating from outside the loading ring radius. Simiu et al. (1984) carried out experiments on 56 small square glass plates in double ring bending. The plates had the nominal dimensions 179 × 179 × 6 mm 3 . The mean thickness was 5.4 mm. The glass was new in the as-received condition and it had been obtained from the same manufacturer and batch. The loading ring radius was 25.4 mm and the segmented reaction ring radius 60.3 mm. All specimens were subjected to ramp loads that generated an average rate of stress of 0.8 MPa/s inside the loading ring. The load-duration until failure ranged from 48 to 117 s. It is not known whether it was the tin side or air side of the glass plates that was subjected to tensile stress. The tin side is defined as the side of the glass that was in contact with the molten tin bath in the float process production method.
This experiment is selected for a number of reasons. The data report is complete with values for the fracture stress even when the failure originated from outside the loading ring. Because Simiu et al. (1984) reported the fracture locations it is possible to make comparisons with the simulated failure origin data. The data is challenging to model. A Weibull distribution for the strength can be rejected, cf. Sect. 6.1. The modelling of the surface flaws condition is simplified when edge failures are eliminated.
Using a formula for a flat circular plate of constant thickness, Simiu et al. (1984) calculated the in-plane MPTS for each fractured specimen. The stress was calculated at the centre of the plate, even when the fracture origin was not located within the loading ring radius. Twelve of the data points, however, were associated with failures originating from outside the loading ring. Those values have been readjusted by this author in order to reflect the MPTS at the actual failure location rather than the MPTS inside the loading ring. The adjustments were made based on finite element calculations with the computer software ABAQUS/CAE (2013). The loading rings were modelled by analytic rigid surfaces. The glass part was modelled with 20node quadratic solid elements with reduced integration, although it would also be possible to use continuum shell elements. The number of through-the-thickness elements was 5 and the number of elements in the plane was about 9500. Only a quarter of the plate was modelled for symmetry reasons. It was assumed that Young's modulus is 70 GPa and Poisson's ratio is 0.23 (Bourhis 2008). A friction coefficient of 0.1 was used in modelling the contact between loading ring and glass parts.

Results
Virtual glass specimens were tested until failure in double ring bending and the results were compared with data from the experiment conducted by Simiu et al. (1984). The analysis was carried out using the software ABAQUS/CAE (2013) and MATLAB (2016). The following cases were investigated, viz. a single population of Pareto distributed flaw sizes using either the MPTS mode I fracture criterion or oblique angle mode I criterion or the mixed mode criterion, and a two-population concept for the flaw sizes using only the MPTS mode I fracture criterion.

Single population of flaws
A fracture stress distribution was simulated based on the oblique angle mode I fracture criterion, cf. Eqs. (1), (2) and (4) while supposing that the surface condition is characterized by a single population of Pareto distributed flaw sizes. This is illustrated in Fig. 3a and the values shown are the MPTS at the failure origins. The histogram in Fig. 3a is normalized so as to reflect a probability density function. The area of each bar is the relative number of observations. The total sum of the bar areas is less than or equal to 1 depending on whether or not some of the data lies outside the bin limits. The sampled distribution was compared with a Weibull distribution and the goodness-offit was tested using the Anderson-Darling (AD) statistic (D'Agostino and Stephens 1986). No significance was obtained in a test at the 5% level. The simulated distribution appears to be indistinguishable from a Weibull distribution. An ordinary Weibull distribution was fit- ted to the empirical data set using the maximum likelihood method and the estimated parameter values were k = 78 MPa and m = 3.8. Using the AD statistic it could be rejected at the 0.4% level that the empirical data set is Weibull distributed. The simulated distribution was optimized so as to match the ordinary Weibull model which was fitted to the empirical data.
The optimization was carried out by varying the underlying Pareto distribution parameters until the simulated strength distribution was similar to the Weibull distribution that was fitted to the empirical data. The similarity was measured by fitting a Weibull distribution to the simulated sample and comparing the so fitted Weibull parameters with the parameter estimates of the Weibull model that was fitted to the empirical data set. See Fig. 3a where both the empirical data set (red bars), the fitted Weibull density function (solid line) and the simulated distribution (black bars) are illustrated. The strength distribution was further simulated using the MPTS mode I criterion and the mixed mode criterion. The Pareto parameters were selected so that the strength distribution could be fitted by a Weibull distribution with scale and shape parameters k = 78 MPa and m = 3.8. Table 1 contains the Pareto parameter values so far discussed according to the three fracture criteria. In Fig. 3b the simulated fracture locations in the radial direction are shown together with the empirical for the oblique angle mode I case. Using a twosample AD test (Scholz and Stephens 1987;Trujillo et al. 2007), a significant deviation between the pair of data sets could be detected. The spatial distribution of failures is further illustrated in Fig. 4a.
The critical stresses in a single cut out specimen are obtained by solving for σ n in Eq. (2) after substituting for the fracture toughness value in Eq. (1) as illustrated in Fig. 4b.
Considering the various fracture criteria, the following was noted while using identical Pareto parameter values for generating the strength distribution. The mode I fracture criterion in the oblique angle case produced only a very small difference in the strength data sample compared with the mixed mode criterion, cf. Eqs. (1), (2), (4) and (5). The 63rd percentiles deviated by less than 3%. However, taking mode II shearing displacement into consideration increased the proportion of failures originating from outside the loading ring by 20%. Comparing the flaw-orientation independent MPTS mode I criterion, Eqs. (1), (2) and (3), with either of the two other criteria yielded a significant difference in the data samples; the 63rd percentile of the simulated strength was more than 10% lower while using the MPTS case. The proportion of failures originating from outside the loading ring increased by over 60%. The results are illustrated in Fig. 5.

Two populations of flaws
It is possible to obtain a simulated distribution like the one shown in Fig. 3c while assuming that the flaws originate from two different populations, see Sect. 3. The flaw model parameters are given in Table 2. The resulting distribution could not be distinguished from the empirical data set with any statistical significance at the 5% level judging from the two-sample AD test statistic ( p = 0.64). In Fig. 3d the simulated fracture locations are shown together with the empirical. It was found using the same test statistic that a significant departure exists from the hypothesis that the experimental and simulated fracture location data sets come from equal distributions.

Discussion
Providing for consistency in a glass failure prediction model calls for its foundation to be laid on physically sound concepts such as the WLP. The WLP captures an essential feature of brittle material failure. The existence of Griffith flaws is another physical concept to build upon. The Weibull distribution implements the WLP which makes it an attractive choice for a glass failure prediction model, at least from a theoretical point of view. All major standards including the European draft of a Eurocode of glass acknowledge Weibull's Eq. (8) in one form or another (prEN 16612:2013). A number of studies, however, have indicated that the Weibull distribution does not provide a superior fit compared with a lognormal or normal distribution (Lü 1997;Calderone et al. 2001;Veer et al. 2009;Huerta et al. 2011;Kinsella and Persson 2016). It has been noted that the estimated value of the Weibull shape parameter varies quite significantly from one sample to another in experiments (Ritter et al. 1985;Carre 1996;Huerta et al. 2011). Some researchers have called for abandoning   (Calderone et al. 2001). But to adopt a normal distribution in this case is to favour with a model lacking in failure-based physical concept. In contrast, by using the numerical method in this paper it is possible to keep intact the WLP as well as the Griffith flaws assumption while producing data fits equal or superior to the Weibull model. Figure 3a illustrates that it is possible to simulate a Weibull distribution using this numerical tool while assuming that the sur-face flaws are sampled from one single population of Pareto distributed sizes. In keeping with recent experimental findings using optical scanning techniques, cf. Wereszczak et al. (2014), it was assumed that the flaw density is 2 cm −2 . From a theoretical point of view, the shape parameter of the simulated Weibull distribution should relate with the Pareto shape parameter according to Eq. (11) if the stress state is uniform uniaxial. At any rate, the Weibull distribution doesn't actually model the experimental data that was compared with. While exploring the possibility of implementing two flaw populations, the idea is to distinguish between one large flaws population of flaws greater in depth than about 100 microns and one smaller flaws population. An idea along a similar line was proposed by Mencik (1992), cf. Sect. 3. The purpose is to model the bimodality that is frequently encountered in the strength distribution from practical experiments. Turning to Fig. 3c it is evident that an acceptable fit can be achieved with a two-population flaws concept. Moreover, this fit is at least as good as the fitted Weibull model in Fig. 3a as can be seen by comparing the p-values from the AD tests. Moreover, the tail of the distribution is important when calculating the design value. Therefore, when choosing between the simulated distributions as seen in Fig. 3a and c, as a matter of fact, the ordinary Weibull distribution appears to provide the most conservative approach. With a two-parameter Weibull distribution, only two parameters are fitted to the data. With the numerical model presented in this paper, the Pareto and Fréchet distributions each require two parameters. As the number of parameters increase, it is only logical that a better fit might be produced. Therefore, the outcome while comparing Fig. 3a with c is rather predictable. However, if it were possible to estimate some of the surface flaw parameters a priori, the numerical modelling tool would gain in potential. Then, these parameter estimates would be based on the material physics. There is a need for more data on the surface flaws condition in glass. Up to date, the published data is scarce. As new techniques become available for examination and assessment of the surface condition in glass, more reliable input data will likely become available for use in this kind of numerical prediction tool. Moreover, with this numerical tool it is possible to simulate the distribution of fracture locations. The simulations were not quite able to model the empirical distribution of fracture location. This is due to the lower mean value in the simulations as well as the longer tail, cf. Fig. 3b and d. However, it may also be due to the fact that a large number of fractures in the empirical data set occured at the loading ring contact circle. About one in five specimens failed under the loading ring. This could have an impact not only on the failure location statistics but also on the fracture stress statistics.
The Weibull shape parameter value that was estimated based on the double ring bending experiment carried out by Simiu et al. (1984), i.e. m = 3.8, might indicate a high dispersion for the experimental data because the value is quite low. The data refers to an experimental campaign carried out almost 40 years ago. In order to investigate the dispersion, a table was organized, cf. Table 3, which contains the estimated Weibull shape parameter values from a range of experiments with the double ring bending device. All listed items in Table 3 refer to experiments on new, annealed glass that was tested in ambient conditions. The experiment carried out by Simiu et al. (1984) is included in Table 3 where the estimated shape parameter value was based on the original data, unadjusted with respect to the true stress at the fracture origin, see also Sect. 5. The results found in Table 3 show clearly that there appears to be nothing unusual about the estimated Weibull shape parameter value in the experimental data of Simiu et al. (1984). However, it is possible that the estimated shape parameter values in Table 3 were affected by the following circumstances. The table comprises both the results from glass tested with the tin side in tension and glass tested with the air side in tension. It has been noted that some experiments with the double ring bending device generate a substantial number of fractures near the loading ring contact area, see e.g. Simiu et al. (1984). Reid (2007) studied the proportion of failures occuring inside and outside the loading ring in coaxial double ring bending tests and compared the observed results with what might be expected based on theoretical considerations using Weibull statistics. He found that a series of 59 small specimens of annealed plates 6 mm in thickness produced anomalous results. The proportion of failures occuring outside the loading ring was substantially greater than expected. Reid hypothesized that this might be related to the glass having to withstand shear stresses outside the loading ring. Due to the equibiaxial state of stress within the loading ring, shear stresses are not present there. Our results show that if a uniformly distributed flaw orientation is considered in the fracture criterion, then there is a significant effect on the observed proportion of failures originating from outside the loading ring while taking mode II shearing displacement into consideration. The proportion increases by 20% with the mixed mode failure criterion. However, disregarding flaw orientation altogether in the fracture criterion, i.e. considering only the MPTS, yields the highest proportion of failures originating from outside the loading ring. Our results are therefore not conclusive with respect to Reid's hypothesis. It depends on whether or not it is assumed that flaw orientation matters. More experiments need to be carried out in order to verify or disprove this hypothesis while taking note of the fracture statistics of failures occuring outside the loading ring.
Although the simulations are more time-consuming than fitting a standard statistical distribution, significant improvements in computational efficiency can certainly be made. There is mounting evidence in the literature, see e.g. Veer (2007), that the fitted models for glass fracture data in general are lacking in potential when using a standard distribution such as the Weibull or Normal distributions. The present study was undertaken in order to explore a novel approach towards the failure prediction of glass. In order to further validate this method, more experiments could be carried out and the surface condition of glass should be investigated further.
The effects of stress corrosion on the strength of glass were neglected in this study. In a future paper, the implementation of subcritical crack growth into the numerical method will be considered.

Conclusions
Using a numerical simulation tool based on the weakestlink principle and assuming the existence of Griffith flaws it is possible to simulate a Weibull distribution for the strength of glass. The incorporation of mode II shearing displacement into the fracture criterion has only a very small impact on the simulated strength distribution when the glass is subjected to double ring bending. In the case of small plates in double ring bending where edge failures can be neglected, it is feasible to model the strength based on a large-flaws and a smallflaws concept while capturing a bimodality in the data set. There is a need for more knowledge and data on the surface condition in glass.