Relationship between strain energy and fracture pattern morphology of thermally tempered glass for the prediction of the 2D macro-scale fragmentation of glass

This work deals with the prediction of glass breakage. A theoretical method based on linear elastic fracture mechanics (LEFM) merged with approaches from stochastic geometry is used to predict the 2D-macro-scale fragmentation of glass. In order to predict the fragmentation of glass the 2D Voronoi tesselation of distributed points based on spatial point processes is used. However, for the distribution of the points influence parameters of the fracture structure are determined. The approach is based on two influencing parameters of fragment size δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} and fracture intensity λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}, which are described in this paper. The Fragment Size Parameter describes the minimum distance between the points and thus the size of a fragment. It is derived from the range of influence of the remaining elastic strain energy in a single fragment taking into account the LEFM based on the energy criterion of Griffith. It considers the extent of the initial elastic strain energy U0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U_0$$\end{document} before fragmentation obtained from the residual stress as well as a ratio of the released energy η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} due to fragmentation. The Fracture Intensity Parameter describes the intensity of the fragment distribution, and thus the empirical reality of a fracture pattern. It can be obtained by statistical evaluation of the fracture pattern. In this work, the fracture intensity is determined from the experimental data of fracture tests. The intensity of a fracture is the quotient of the number of fragments in an observation field and its area and is assumed to be constant in the observation filed. The fracture intensity and the correlation between a constant intensity and the Fragment Size Parameter was determined. The presented methodology can also generally be used for the prediction of fracture patterns in brittle materials using a Voronoi tesselation over random fields.

tribution, and thus the empirical reality of a fracture pattern. It can be obtained by statistical evaluation of the fracture pattern. In this work, the fracture intensity is determined from the experimental data of fracture tests. The intensity of a fracture is the quotient of the number of fragments in an observation field and its area and is assumed to be constant in the observation filed. The fracture intensity and the correlation between a constant intensity and the Fragment Size Parameter was determined. The presented methodology can also generally be used for the prediction of fracture patterns in brittle materials using a Voronoi tesselation over random fields.
Keywords Fragmentation · Tempered glass · Fragment size · Fracture intensity · Elastic strain energy · Spatial point process · Voronoi tessellation · Fracture pattern Critical energy release rate K I c Critical stress intensity factor η Energy relaxation factor A η Energy relaxation zone of a fragment U η Elastic strain energy in the relaxation zone A η U 0 Initial strain energy U 1 Remaining strain energy U R Relative remaining strain energy N D Fragment density in the observation field with the length of D

Introduction
Glass is one of the most popular building materials today. However, the tensile strength is governed by small flaws in the surface which reduce the actual engineering strength of annealed float glass to 30-100 MPa ). Due to the residual stress state thermally tempered glass shows a greater resistance to external loads and in case of failure it is quite safe in terms of cutting and stitching due to the small, blurred fragments. Therefore, thermally tempered glass is also known as tempered safety glass. The residual stress state is obtained by the tempering process and is approximately parabolic distributed along the thickness with the compressive stress on both surfaces and an internal tensile stress in the mid-plane. By imposing a compressive residual stress at the surface, the surface flaws will be in a permanent state of compression which has to be exceeded by externally imposed stress before failure can occur (Schneider 2001;Nielsen et al. 2010;Pourmoghaddam et al. 2016;Pourmoghaddam and Schneider 2018a). The amount of the residual surface compressive stress largely depends on the cooling rate and therefore on the heat transfer coefficient between the glass and the cooling medium (Gardon 1965;Aronen and Karvinen 2017). However, thermally tempered glass will fragmentize completely into many pieces, if the equilibrated residual stress state within the glass plate is disturbed sufficiently and if the elastic strain energy in the glass is large enough, Fig. 1. The so-called Rupert's drop with bulbous head and thin tail can withstand high impact or pressure applied to the head, but explodes immediately into small particles when the tail is broken (Silverman et al. 2012). The fragmentation is the direct consequence of the elastic strain energy that is stored inside the material due to the residual stress state. The fragment size depends on the amount of the stored energy. Small fragments are caused by a high stored strain energy due to the high residual stress state originating from the extremely rapid cooling. And lower residual stress states result in larger fragments due to lower stored strain energy (see Fig. 2). Thus, not only the stress but also the thickness of the glass plate plays a role in determining the strain energy.
Thus, for the same residual stress distribution a higher stored energy is determined in the thicker plates. This has been investigated and proved experimentally by several authors (Akeyoshi and Kanai 1965;Barsom 1968;Gulati 1997;Lee et al. 2012;Mognato et al. 2017;Pourmoghaddam and Schneider 2018b).
The fragmentation pattern, i.e. the fracture structure, the fragment size and thus the fragment density are the direct consequence of the amount of the initial strain energy U 0 which is available before the fragmentation sets in and can be written as: where it should be noted that in this study, it is assumed that the glass is only subjected to residual stress from the thermal tempering U 0 = U Residual stress . When the glass is fragmentized the internal energy converts into a number of forms of energies (Reich et al. 2013;Nielsen 2017), e.g. for creating new fracture surfaces (cracking and crack branching), kinetics and sound energy. Hence, the remaining strain energy in a fragment U 1 after the fragmentation sets in must be expressed as: where U other represents other energy consuming effects such as heat. In this study, only the relative remaining strain energy U R as the ratio between the remaining strain energy and the initial energy is considered. The relative remaining strain energy can be written as: The other forms of energy are not included here. The remaining stress state and the resulting remaining elastic strain energy in a single fragment has been calculated numerically by Nielsen (2017). Further experimental and numerical investigations of fragment's deformation and strain energy were carried out by Nielsen and Bjarrum (2017). The change in strain was determined by comparing the surface shape of a fragment before and after fracture.
Several models for relating the fragment size to the residual stress state have been suggested in the literature, e.g. Acloque (1956), Barsom (1968), Gulati (1997), Shutov et al. (1998), Warren (2001) and Tandon and Glass (2005). Some of these works have proposed models for the fragments size based on an energy approach. Most of the works try to establish an analytical model for the fragment size considering the release of the so-called tensile strain energy defined as the part of the strain energy resulting from the mid-plane tensile stress alone. A theoretical method based on fracture mechanics was developed by e.g. Warren (2001), Molnár et al. (2016). The Voronoi tesselation of randomly distributed points in the glass plane leading to a cell structure similar to glass breakage was mentioned in Molnár et al. (2016). However, the statistical distribution method of the points in the plane itself is significant when it comes to the prediction of the fracture pattern. With regard to the parameters of the fracture structure, the comparison of different spatial point processes is therefore important in order to find the best method of point distribution. The motivation is that Voronoi tessellation of points distributed in the plane results in a fracture structure based on both fracture mechanical and stochastic parameters. Hence, the main objective is to determine the significant influence parameters of the glass fracture in order to develop a method for the prediction of 2D macro-scale fragmentation of glass. In this work, two influence parameters of the fracture structure, the Fragment Size Parameter δ and the Fracture Intensity Parameter λ are described and determined. Using these two parameters the spatial point process can be influenced and calibrated by the energy density and empirical reality of a fracture pattern.
In order to determine the Fragment Size Parameter δ a method was developed with which the minimum distance between two points in a point cloud can be determined by selecting the energetic criterion (Griffith 1920) from the linear fracture mechanics with respect to the strain energy state before and after the fragmentation. In other words; the present work aims at determining the minimum distance between randomly distributed points in a plane representing fragments based on the strain energy state which is stored in the glass plate due to the thermal tempering. Subsequently, the fragment density is estimated using Hexagonal Close Packed (HCP) and a uniformly distributed point group.
The Fracture Intensity Parameter λ describes the degree of intensity of the fracture in a limited observation field and is a value for the fragment number within the observation field. This parameter was determined from the experimental data of fracture tests.

Energy conditions
When external forces or residual stresses deform an elastic body, these stresses perform work as their points of application are displaced. This work is stored in the body as elastic strain energy. The total strain energy U stored in a deformed linear elastic, isotropic body is obtained by integrating the energy per unit volume over the volume of the body: where σ i j is the stress tensor and ε i j is the strain tensor. The residual stress in thermally tempered glass is distributed parabolically along the glass thickness t (Fig. 3). This parabolic stress distribution σ (z) can be written in terms of the surface stress σ s as: using symbols defined in Fig. 3. The parabolic stress distribution is in equilibrium and symmetric about the mid-plane. The magnitude of the surface stress is approximately twice the tensile stress (2σ m = −σ s ). The zero stress level is at a depth of approximately 20% of the thickness t. The stress state σ is assumed to be planar. Using cylinder coordinates r for the radius and θ for the polar angle the constitutive relationship for plane stress is defined by: with the Young's modulus E and the Poisson's ratio ν. The equilibrated stress state is also assumed to be hydrostatic for field stresses. A planar hydrostatic stress state means that no shear stresses occur (τ r θ = 0) and the normal stresses are always principal stresses, which are equal in the plate plane (σ r = σ θ = σ (z)) and zero in the direction of the thickness (σ z = 0). Hence, the total initial elastic strain energy U 0 can be written as: Inserting the residual stress field from Eq. (5) and integrating over the cylinder we find: which is the initially stored strain energy in a cylindrical body of the radius R, thickness t and the residual surface stress of σ s . The Eq. (9) can be written in terms of the strain energy per unit surface area of any given base shape of a body by dividing with the base area for the cylinder (Barsom 1968;Gulati 1997;Warren 2001;Nielsen 2017): The energy can also be written in terms of the energy density U D : which is the amount of elastic strain energy stored in the system per unit volume and thus only depends on the residual stress and the material properties.

Voronoi tesselation over spatial point processes
This section provides a brief introduction to the two mathematical tools that play a central role in the latter algorithm: Spatial point processes and Voronoi Diagrams. Further details and exhaustive presentations of these topics are covered in Baddeley et al. (2016), Wiegand and Moloney (2014), Schmidt (2015), Møller (1994), Okabe et al. (1992) and Ohser and Schladitz (2009).

Spatial point processes
A point process describes random configurations of points P = {p 1 , . . . , p n } in a continuous bounded set K (which is in the case of this paper the 2D domain of the glass fracture pattern with area content |K |).
The number of points n is itself a random variable N that typically follows a discrete Poisson distribution. The least complex point patterns exhibit complete spatial randomness (CSR), i.e. the point locations p i occur independently and uniformly (equal likelihood) over the domain K . The probability density of a point pattern P with point process parameters θ and distribution of the point locations f θ n (p 1 , . . . , p n ) is given by Here, P is one of the possible point patterns and f (P) can be interpreted as 'the relative chance of obtaining the point pattern P' (Baddeley et al. 2016). Note, that in Eq. (12), the probability density of a CSR does not depend on the locations of the points p i , i.e. all possible patterns of n points are equally likely, the different points are furthermore independent and each point is uniformly distributed over the domain K . Further point processes can be defined by probability densities f P (P|θ) which differ in their mathematical form from Eq. (12) (i.e. a CSR process), this will serve as the basis for two more point processes presented later in this section.
The behavior of spatial point processes can be characterized in terms of first-order and second-order properties. The first-order property is described by the intensity λ(p) (mean number of events per unit area at the point p) and the second-order property reflects the spatial dependence in the process (clustering or repulsion). More complex spatial point processes are able to enforce either clustering or repulsion of the points, which is observable in the second-order statistics of these processes. Figure 4 shows examples of point patterns which exhibit CSR and two examples violating CSR due to regularity (repulsion) and clustering. Based on inspection of typical fracture patterns, such as shown in Figs. 1 and 2 a repulsion behaviour of the underlying spatial point process can be concluded. Repulsion processes possess a repulsion domain, which is located around each seed point to avoid points being too close to each other. When dim K = 2, this domain is a disk with radius δ HC , which is a model parameter of the underlying spatial point process. This kind of pairwise interaction of the points can be modelled via repulsion/inhibition processes. In the context of this paper Matérn Hardcore Processes (MHC) (Matérn 1960) as well as Strauss (repulsion) processes (SR) (Strauss 1975) are investigated, which belong to the family of the Gibbs processes (Baddeley et al. 2016).
In contrast to the CSR process with one parameter (intensity λ), the HC process possesses two parameters, the intensity λ HC and the (minimal) hardcore distance δ HC , thus it is a two-parametric model with probability density (Baddeley et al. 2016): where c HC is the normalization constant. Matérn (1960) proposed three schemes, called Matérn type-I, type-II, and type-III hardcore point processes, for constructing repulsive point processes. A type-I process is gathered by sampling a primary process from a homogeneous Poisson process with intensity λ, and then deleting all points with distance less than δ HC . The Matérn type-II process assigning each point an 'age', which defines, in case of point-distances less than δ HC , that the older point is kept while the younger point is deleted in the sampling process. A Matérn type-III process lets a newer event be thinned only if it falls within radius δ HC of an older event that was not thinned before. According to Baddeley et al. (2016), the hard core process is a CSR conditional on the event, that it satisfies the hard core constraint. For the modelling of the fracture pattern, this means, that the intensity λ HC is equivalent to the intensity of a CSR λ while only keeping points, that are not closer than δ HC . In the latter of this paper it is highlighted, that the hardcore constraint introduces significant differences in the second-order statistics of a point process and it is further shown, that fracture patterns from experiments are in favour of a hardcore constraint. Acc. to Illian et al. (2008), the intensity of a HC λ HC can be estimated from the intensity of a CSR λ via where V = π · δ 2 HC . Analogously to the derivation of Eq. (29), combination of Eqs. (27) and (14).
In general, the hard core process is an appropriate model given, that it is physically impossible for two points to lie closer than a distance δ HC . If close pairs of points are not impossible but unlikely to occur, the Strauss (repuslion) process (SR) 'penalises' rather than 'forbids' close pairs of points. In the context of this paper Strauss (repulsion) processes SR (Strauss 1975) are investigated, which belong to the family of the Gibbs processes. Without going in detail here, c.f. (Stoyan et al. 1995), the density of a Strauss point process reads where c S R is a constant, β S R > 0 is the 'fertility' parameter, 0 ≤ γ S R ≤ 1 is the interaction parameter, and s(p) is the number of point pairs, of which the Euclidean distance is less than δ HC . The Matérn Hardcore Process is a special case of the Strauss process with γ S R ≡ 0, where γ S R can be interpreted as the inhibition parameter. The considerations before lead to a total of two respectively three parameters, which are necessary to characterise either the Hardcore Process θ HC = {λ HC ; δ HC } or the Strauss Process: θ S R = {β S R ; δ HC ; γ S R }. However, in Sect. 4 it is derived, that the parameters θ HC as well as θ S R are motivated from and are connected to distinct parameters from Linear Elastic Fracture Mechanics (LEFM) considerations. According to Baddeley et al. (2016), β is not the 'intensity' of the HC or SR model, instead β should be seen as the spatially varying 'fertility' which is counterbalanced by the 'competitive' effect of the hard core to give the final intensity, thus the 'intensity' is the product of fertility and competition. β is referred to as the (first-order) trend which is of more interest than the resulting 'intensity' of a point process.
The observed number of seeds in the inspection window n can be interpreted as a further parameter of the processes, however conditionally on n, the processes possess the number of parameters as stated before. The sampling of a point process is in general not a simple task (as e.g. the normalizing constant Z may be unknown as in the Strauss process case), thus usually Monte Carlo methods have to be applied as they do not require normalized probabilities for sampling (Baddeley et al. 2016;Møller et al. 1999).

Voronoi tesselation
For a given configuration of points P in K , called the seeds, the Voronoi cell associated to the seed p i ∈ P, denoted as V (p i ), corresponds to the region in which the points are closer to p i than to any other seed in P: The Voronoi diagram generated by P is the set of the Voronoi cells {V (p 1 ), . . . , V (p n )}. In Fig. 5, the Voronoi tessellation of randomly distributed points is shown. Voronoi diagrams entirely partition the domain K without region overlap. By using the Euclidean distance (q = 2) in Eq. (16), Voronoi cells are guaranteed to be convex polygons. Refer to Okabe et al. (1992) for a thorough treatment of Voronoi tessellations and their properties.
Matérn Hardcore and Strauss repulsion processes in the context of this paper are used to generate the The basic idea for the theoretical prediction and simulation of the fracture pattern is to distribute points with a certain distance δ to each other based on the locally acting stress respectively the stored elastic strain energy, and with a fracture intensity of λ obtained by evaluation of fracture patterns in the plane. The final fracture pattern is created by Voronoi tessellation of the resulting point cloud. Points can be distributed using different distribution methods (see Sect. 3). Taking into account the size and intensity parameter, the following point distribution dependency results as the estimation parameters for the theoretical prediction of the 2D macro-scale fracture structure: This work focus on the elastic strain energy resulting from the residual stress state from the thermal tempering process. In order to predict the fragmentation and for the recognition and the generation of the fracture The Fragment Size Parameter δ is equal to the hardcore distance of the hardcore and Strauss spatial point process and thus can be understood as a local minimum point distance. For the determination of the Fragment Size Parameter δ not only in dependence of the residual stress based on the thermal tempering, but on any kind of stress (e.g. stress resulting from an external load), the linear elastic fracture mechanics (LEFM) based on the energy criterion introduced by Griffith (1920) is used. Part of the energy is released by new surfaces generating from cracking and branching of progressive cracks. The fracture pattern can be estimated with the help of the energy concept in fracture mechanics, taking into account stochastic fracture pattern analyses. This method is compared with the results of the finite element simulations in Nielsen (2017).
In high-speed images (Nielsen et al. 2009), it was observed that so-called "whirl-fragments" were generated by a whirl-like crack propagation. It was also observed that progressing cracks branched at an angle of 60 and formed a hexagonal fracture. A hexagonal fracture structure is assumed to be the "perfect" fracture pattern or "perfectly" broken glass plate. Predicting such a "perfect" fracture structure is easy, provided δ is known. All points in the plane have the same distance to each other and can be distributed by Hexagonal close packing (HCP) of points (see Fig. 6a). The Voronoi tessellation takes place via the Delaunay triangulation of the HCP-distributed points (Fig. 6b) and results in the cell structure of a honeycomb (Fig. 6c). Thus the number of fragments in an observation field can be predicted well.
However, since the fracture pattern will not be perfect in reality (see Fig. 7), due to stochastic material properties and stress distribution from the production process, another parameter λ, which can be obtained from experimental data and describes the intensity of the fragment distribution, and thus the empirical reality of a fracture pattern, must be added to the method.
The Fracture Intensity Parameter λ, which is a characteristic value for the number of fragments in an observation field, was determined using the experimental data of the fracture tests in carried out in Pourmoghaddam and Schneider (2018b). In Fig. 7, three samples of fracture patterns are shown for the same glass thickness but different residual stresses respectively elastic strain energy density U D . Fracture intensity is determined by placing 50 mm × 50 mm observation fields on the fracture patterns and determining the average number of fragments in correlation to the energy density for each sample (Pourmoghaddam and Schneider 2018b).

Fragment Size Parameter δ
The determination of the Fragment Size Parameter δ is based on the energy conditions described in Sect. 2. As it is shown in Fig. 6a, the Fragment Size Parameter δ = 2r 0 is the distance between two neighboured points. For the development of the approach a HCP distributed point cloud is assumed. This is necessary because we do not want to go into the intensity of the fracture pattern at first and want to describe the point distance purely physically. There are three assumptions which were necessary for the determination of the Fragment Size Parameter respectively for the minimum distance of two distributed points.

Assumption 1
The first assumption is that the glass plate will break into cylindrical fragments. It was assumed that the sphere of influence for the stored elastic strain energy U of each point representing a fragment is a circle, which is here called "Energy circle" respectively a cylinder in 3D with a radius r 0 before fragmentation. As can be seen in Fig. 6a, the energy circles touch but do not overlap. Consequently, all energy in the observed field is distributed in these energy circles. Thus, each radius r 0 depends on the elastic strain energy in the influence area of the respective point. The free gaps between the energy circles are zero energy areas and have no influence on the further calculation. This is legitimate, as we are distributing the total strain energy in the observed field through the energy circles.
Assumption 2 In the case of fracture the stored elastic strain energy will release and there will be a relaxation of energy in the fragment. However, we assume that in the case of fracture the elastic strain energy is not relaxed completely but just by a relaxation factor η: This means, that there will be a remaining strain energy U 1 in the fragment after the fragmentation. This was shown numerically in Nielsen (2017) by FE simulations on a fragment before and after the fracture sets in. For the energy sphere described under Assump- Fig. 8 Energy circle before fragmentation with radius r 0 and after fragmentation with radius r 1 ; A η is the surface of the energy relaxation zone tion 1, this means that the energy circle shrinks to a smaller circle with a radius r 1 after fragmentation (see Fig. 8). The area in which the energy relaxes is called the Energy relaxation zone. The area of the Energy relaxation zone A η can be easily calculated as follows: In addition, the elastic strain energy U η in the relaxation zone can be calculated as follows: whereinσ describes the stress state as a factor of the integrated stress function through the thickness.σ can be determined asσ m = 4 5 σ 2 m by calculating with the residual tensile mid-plane stress andσ s = 1 5 σ 2 s for the residual compressive surface stress for a parabolic stress distribution from tempering. U η is therefore the elastic strain energy released during fragmentation.
Assumption 3 It is assumed that in the event of a fracture, the elastic strain energy is released through the generation of new fracture surfaces. Other forms of energy such as acoustic or heat are neglected. Therefore, we assume that all released elastic strain energy U η is converted to the fracture surface energy Γ : where A f r = 2πr 1 tρ * is the fracture surface by assuming a cylindrical fragment. ρ * is a correction factor of the fracture surface, since the fracture surface is not flat in reality but irregularly shaped (Quinn 2016).
Here A f r,act is the actual fracture surface and A f r,cyl is the lateral surface of a cylinder. Therefore, for an assumption of a cylindrical fragment ρ * = 1. γ is the specific fracture surface energy and can be described in terms of the critical energy release rate: Depending on the stress state, we obtain: , Plane strain (24) where in K I c is the critical stress intensity factor. Inserting the Eq. (20) in Eq. (21) and resolving it to r 0 , we obtain: Furthermore, the relative remaining elastic strain energy U R,Rem can be described as: According to the first assumption, the glass plate will break into cylindrical fragments. In Fig. 9 the relative remaining elastic strain energy U R,Rem in correlation with the area of the base shape of the fragment A f r,base normalized by t 2 is shown in comparison with the FE results in Nielsen (2017). As can be observed in Fig. 9, the analytical results of the fragment size from a value of approximately A f r,base /t 2 ≥ 1.5 corresponds to the FE results. The difference in the case of very small fracture surfaces is due to the fact that the FE simulation was carried out on a fragment volume and therefore a multi-axial stress state was present. In the FE simulations stresses in the thickness direction occurred in the fracture state, whereas the analytical solution was based on a continuous plane stress state. Now we have described the distance between the distributed points δ = 2r 0 , which directly affects the fragment size and so the fragmentation density N , in dependence on the residual stress respectively the elastic strain energy U . The method significantly depends on the relaxation factor η (see Fig. 10), which can have values from 0 to 1. The relaxation factor of η = 0 would mean that there is no fracture and η = 1 that the energy is relaxed completely. Thus the fragment size or the minimum distance between two points for a hardcore spatial point process is mainly dependent on the energy relaxed in the fracture state.
As described before and shown in Fig. 10, the fragment density is affected by the relaxation factor η. The greater η, the higher the percentage of the stored elastic strain energy that is released during the fragmentation. A larger energy produces more crack surfaces and thus also a finer fracture pattern or in other words a larger fragment number within an observation field. For example for a residual surface compressive stress of 100 MPa we calculate in an observation field of size 50 mm × 50 mm a fragment density of N 50 = 14 for a relaxation factor of η = 0.05 and N 50 = 132543 for a relaxation factor of η = 0.9.
The Fragment Size Parameter δ was determined using the experimental data of fracture tests. In Pourmoghaddam and Schneider (2018b) specimens with different thicknesses were thermally tempered with different degree of tempering so that specimens with different residual stresses were obtained for the fracture tests. After the fracture tests the average fragment weight was determined from more than 130 fragments per specimen. Knowing the density of glass (2500 kg/m 3 ) and the glass thickness t the area of the base shape of the fragment (fragment volume divided by thickness) can be recalculated from the weight. Using the energy density U D for the determination of the correlation, the curves of base area and thus the Fragment Size Parameter δ for the different thicknesses coincide to one line. In Fig. 11, the elastic strain energy density calculated from the different residual stresses of the specimens is applied over the Fragment Size Parameter δ. Each of the black triangles represents the average of more than 130 fragments per specimen. For the calculation of δ from the experimental data a fragment with the number of edges n → ∞ (cylindrical fragment) is assumed which contains other edge numbers and the Fragment Size Parameter δ is determined over the radius of the fragment. The red line represents the fragment size calculated analytically using Eq. (25). For the analytical calculations the measured residual stress values of the specimens used in Pourmoghaddam and Schneider (2018b), a Young's modulus of E = 70000 MPa, a correction factor of ρ * = 1 and the plane stress state have been taken into account.
As it can be observed in Fig. 11 the analytical curve fits for a relaxation factor of η = 0.123. A relaxation factor of η = 0.123 means that the relative remaining elastic strain energy U R,Rem = 0.77 and thus approximately 77% of the initial elastic strain energy U 0 remains in the fragment and correspondingly 23% of the initial energy releases during the fracture process. In Fig. 12 the correlation between the elastic strain energy density and the fragment density for the initial energy density U D,0 , remaining energy density U D,1 and the released energy density U D,η . Thus a direct correlation between the energy in pre-fracture, fracture and post-fracture state is shown.

Deterministic fragmentation process
The fracture intensity is another important parameter for the theoretical prediction of the fragmentation of glass. The fracture intensity is a characteristic value for the fragmentation behavior contains information about the fragment density (Number of fragments in a given area) of the fracture structure. From a deterministic Fracture Mechanics point of view, the material can be expected to be homogeneous and thus the spatial location of fragment centers does not depend on the actual location within an object under investigation. As a starting point for an estimation of the locations of the centers of the glass fragments this assumption serves well and is called 'deterministic fragmentation process'.
For the computation of the fragment density N D in an observation field a hexagonal fracture structure with a fragment side length a = 2r 0 / √ 3 is now assumed (Fig. 6c). In a square observation field with the side length D the fragment density N D can be described as: Under the assumption homogeneity of the point process the intensity λ is constant over K which further implies, that the expected number of points falling in an observation region with the side length D is proportional to its area: In order to obtain hexagonal Voronoi cell structure (honey comb), points can be distributed by Hexagonal Close Packing (HCP). Combining Eqs. (27) with (28), an approximation for the expected fracture intensity λ HC P of a honey comb, motivated from deterministic fracture mechanics, can be expressed in terms of the Fragment Size Parameter δ = 2r 0 as: In Fig. 13a, a double logarithmic interrelation between fracture intensity of the honey comb and fragment size is shown. The fracture intensity decreases with a coarser fracture structure.λ HC P can further be rewritten in terms of the energy relaxation factor η as: In Fig. 13b, a double logarithmic interrelation between fracture intensity of a honey comb and the energy relaxation factor is shown. The fracture intensity increases with higher energy relaxation.
In Fig. 14a-c the Voronoi tesselation of HCP distributed points with the point distance δ is shown for a plate with the residual mid-plane tensile stress of 50 MPa. The respective fragmentation density is shown  Fig. 14a to η = 0.2 in Fig. 14c. It is shown that for the same residual stress of 50 MPa the point distance decreases with higher values for η.
The assumption of a hexagonal fracture structure is eased to the assumption of a polygonal fracture structure induced by the Voronoi Tesselation over a Hardcore Process (HC). The fragment density N D is now connected to the spatial point process intensity. Analogously to Eq. (29), using Eq. (14) and combining it with Eq. (28), an improved estimation of the expected Fracture Intensity Parameterλ HC can be derived in terms of the Fragment Size Parameter δ = 2r 0 as: However, despite Eq. (32) takes into account the magnification of the observed intensity from the simulations to estimate the underlying HC process intensity, the connection to the HCP is still basis of the derivation of the estimator.

Stochastic fragmentation process
If fracture tests are conducted, it can be observed, that the center locations of fragments do vary stochastically within an object under investigation due to different reasons such as inhomogeneous thermal pre-stressing of a glass pane, inhomogeneous load application etc. Analogously to Sect. 4.3.1, the same investigation was carried out but this time using random point locations in an observation field of size 50 mm × 50 mm with a uniform distribution for the point locations instead of HCP. This kind of randomly distribution led to a varying distance between points. However, the distance determined on the basis of the relaxation factor η was given as the minimum distance for the point process.
In Fig. 15a-c the Voronoi tesselation of uniformly distributed points with the minimum point distance δ min is shown for a plate with the residual mid-plane tensile stress of 50 MPa. In comparison to the fragment density for the case of HCP distributed points, the number of fragments is underestimated for uniformly distributed points. Until η = 0.04 the difference is not significant yet. However, the difference is more pronounced for the higher relaxation factor of η = 0.2. This is because on the one hand, Eqs. (29)-(31) are derived as approximations for the underlying stochastic point process as prior expectations motivated from the deterministic fracture process and on the other hand, in the deduction of Eqs. (29)-(31) a minimum energy requirement was not introduced for the randomly distributed points.

Experimental determination of the Fracture Intensity Parameter λ
The Fracture Intensity Parameter λ has been determined experimentally from fracture tests on thermally tempered glass specimens. In Pourmoghaddam and Schneider (2018b), 72 glass plates of size 1100 mm × 360 mm and three different thicknesses of 4 mm, 8 mm Intensity Parameter λ (1/mm 2 ) determined from fracture tests (Pourmoghaddam and Schneider 2018b) and 12 mm were thermally tempered with different heat treatment (Pourmoghaddam and Schneider 2018c) for achieving different residual stresses and then fractured according to EN 12150-1. Subsequently, the average fragment density N D of eight square observation fields with a side length D = 50 mm taking a boundary of fragmentation analysis into account, as shown in Fig. 16, was determined by counting the fragments of each observation field. The fragmentation was influenced by the impact in an area under the impact position (impact influence zone). Therefore, only the fragments in the areas under an assumed angle of 45 • left and right from the impact point were considered for the investigations. Now assuming a constant intensity in each observation field the average Fracture Intensity Parameter λ can be calculated for each specimen using Eq. (28).
Using the energy density U D for the determination of the correlation, the curves of the fragment density N 50 and thus the fracture intensity λ for the different thicknesses coincide to one line. In Fig. 17, the correlation between the elastic strain energy density U D , calculated from the measured residual stresses of each specimen, and the Fracture Intensity Parameter λ is presented. It can be observed that the accuracy of the experimental results of the fracture intensity decreases with lower energy density respectively for larger fragments.

Conclusion
Based on fracture tests on glass plates with different heat treatment for achieving different residual stress levels (Pourmoghaddam and Schneider 2018c) this paper deals with the determination of Fracture Structure Parameters for the prediction of the 2D macroscale fragmentation of glass. The Voronoi tesselation of distributed points using calibrated spatial point processes can be used in order to simulate the glass fracture structure. The probability density function of the Hardcore Process, e.g., possesses two parameters, the intensity of the points in an observation region and the (minimal) hardcore distance between point pairs, thus it is a two-parametric model. In this work, due to comprehensive investigations of fracture patterns and analytical considerations based on the LEFM the parameters characteristic for the fracture intensity called Fracture Intensity Parameter λ and the hardcore distance between point pairs called Fragment Size Parameter δ have been determined. These two parameters are applied over the elastic strain energy U D in Fig. 18. The fracture mechanical parameters were laid and the connections to the parameters of the spatial point processes were highlighted. Under assumption of a cylindrical fragment with the radius r 0 of the base area, a Fragment Size Parameter δ was generated as a function of the elastic strain energy remaining in the fragment after the fragmentation considering the relaxed energy. In this work, the energy relaxation factor η was fitted to the results of the fracture tests. The minimum distance between two points (Hardcore distance δ = 2r 0 ) can be determined on the basis of the LEFM. The number of fragments or the number of points representing fragments is taken into account in an observation field K as a function of the temper stresses. The Fracture Intensity Parameter λ to consider the intensity of the fracture structure was also discussed. A correlation between the elastic strain energy density and the intensity of the fracture pattern was established based on the experimental data of fracture tests. An estimator for the constant Fracture Intensity Parameter was derived for two cases: the deterministic Hexagonal Close Packed (HCP) setting as well as a Hardcore stochastic point process, both acting on an observation field K . The generation of fracture patterns under the assumption of HCP leads on the one hand to a visually observable poorer depiction of the glass fracture pattern and on the other side a poorer estimate of the fracture intensity compared to the Voronoi tesselation of the observation domain over spatial randomly distributed point locations p. Thus, further research will head in the direction of assuming point locations for the Voronoi tiling as random variables. The underlying probability distribution has to be estimated from fracture experiments.

Outlook and future research
In order to predict the 2D macro-scale fragmentation of glass considerations from Linear Elastic Fracture Mechanics (LEFM) with approaches from stochastic geometry in terms of spatial point patterns have to be merged. This method is called "Bayesian Reconstruction and Prediction of Glass Fracture Patterns (BREAK)". In Pourmoghaddam et al. (2018a, b), numerical studies on the probability distributions of for different geometrical properties of tesselation over the CSR, MHC and SP with different parameter settings were conducted and evaluated. As no analytical expressions for the posterior distribution of e.g. the number of edges, the area content or the circumference of a typical Voronoi cell of a MHC or SP exist and these quantities are of interest in a fracture mechanical setting, numerically intensive Monte Carlo simulations of these processes were used to numerically infer the respective distributions to allow further fracture mechanical processing of these information. Besides the fracture mechanical relevance a general mathematical value can be addressed to the obtained results, as the posterior distributions for the geometrical properties can be used in other applications without loss of generality.
An example for the estimation of the spatial point process intensity and thus fragment density is depicted in Fig. 19 (note, that in this computation no edge correction was applied as the choice of a correction method itself is non-trivial, c.f. Baddeley et al. (2016);Møller (1994). In Fig. 19 the estimation of the intensity λ HC of the underlying point process as well as the induced Voronoi Tesselation for the fracture pattern of one of the tested glass panes is shown as 3D view and a top view. The process calibration as well as simulation of further fracture patterns will be highlighted in part two of this paper in depth. In order to determine the statistical values of the fracture structure such as fragment edge number, the fragment perimeter, fragment base area, etc. first the fracture image has to be recorded and morphologically processed. Then a spatial point process model fed with the informations from the fracture mechanics (fragment size respectively the hardcore distance between point pairs) and the experimental data of fracture tests (fracture intensity) has to be matched to the fracture image and calibrated to evaluate a candidate model. Subsequently, the fracture pattern can be simulated. The overall methodology with the incorporated theories is depicted in Fig. 20. The results presented in this paper, the incorporated theories as well as further experimental investigations of the fracture structure will flow into future research in the field of 2D macroscale fracture structure prediction. An upcoming pub-  lication deals with the further deduction, implementation and calibration of the method "BREAK" in order to simulate fracture patterns which conserve different statistical properties of the underlying glass fracture patterns such as the distribution of interpoint-distances or the distribution of remaining fracture particle area. The assessment of different spatial point processes as well as the properties of the induced Voronoi tesselations have to be studied and compared against the empirical statistics obtained from the images of the fracture tests on glass plates. The formalisation of the findings in terms of parameter tables as basis for the simulation of the spatial point processes with their Voronoi tesselation will be provided in the second part of this paper.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.