# Abraded Glass Strength: An Ad Hoc Fitting Protocol Based on the Change of Variable Theorem

Authors: Camilla Ferretti; Piero Ganugi; Gabriele Pisano; Francesco Zammori
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This work tackles the problem of finding a suitable statistical model to describe relevant glass properties, such as the strength under tensile stress. As known, glass is a brittle material, whose strength is strictly related to the presence of microcracks on its surface. The main issue is that the number of cracks, their size, and orientation are of random nature, and they may even change over time, due to abrasion phenomena. Consequently, glass strength should be statistically treated, but unfortunately none of the known probability distributions properly fit experimental data, when measured on abraded and/or aged glass panes. Owing to these issues, this paper proposes an innovative method to analyze the statistical properties of glass. The method takes advantage of the change of variable theorem and uses an ad-hoc transforming function to properly account for the distortion, on the original probability distribution of the glass strength, induced by the abrasion process. The adopted transforming function is based on micromechanical theory, and it provides an optimal fit of the experimental data.


---

## Body

## 1. Introduction

The increasing demand for transparency in architecture has led to a step change in the structural role of glass, whose applications have moved from simple windows to load-bearing components. This poses serious engineering problems, since the mechanical response of glass, especially its tensile strength, is governed by the presence of flaws on its surfaces. Unfortunately, such flaws are unavoidable, as they are a direct consequence of both the manufacturing and of the subsequent cutting process. Also, perhaps more importantly, the defective scenario can change over time, due to the action of many abrasion and/or corrosion phenomena, which can increase and alter the number and shape of the cracks on the external surfaces. Consequently, the probability distribution function (pdf) of the tensile strength of a glass-made structural element is not fixed but evolves during its lifetime and may significantly deviate from the pdf of the original pristine material (glass as delivered after quality control).The main issue is that, due to the random nature of shape, size, and distribution of surface flaws, the structural strength cannot be considered deterministic and must be treated in a statistical way [1].In recent years, the mechanical characterization of glass has attracted much academical interest. However, most researchers focused on pristine glass, whereas just a few of them tackle the problem of finding a suitable statistical model for aged glass. This is a major gap of technical literature, because structural elements should assure a very low failure probability (during their whole lifetime), as defined by the Eurocode 0 [2]. Many international standards and guidelines provide corrective factors (or coefficients) that structural engineers can apply to properly account for the effects of aging. However, the reliability of these coefficients is rather questionable as they were empirically defined based on past field experience, even if the use of glass in components with an important structural role is very recent, and the know-how about how glass strength may deteriorate with aging is still very limited. Hence, there is a need for a refined statistical model to recalibrate the corrective coefficients (accounting for aging) in structural design.In the literature, glass strength is generally modeled using a 2-parameter Weibull probability distribution [3, 4]. This distribution, in fact, fully complies with the mechanical behavior of glass beams or plates, whose macroscopic collapse is always due to the unbounded growth of only one crack (weakest-link-in-the-chain concept). Nonetheless, such distribution has been recently questioned due to its inability to fit experimental data at the lower quantiles [5], and a better fit has been obtained by substituting the 2-parameter Weibull with a 3-parameter Weibull or with a left-truncated Weibull (LTW) distribution [6]. In both cases, it is assumed that collapse cannot occur below a minimum threshold limit of the tensile stress, and so a lower limit is introduced in its pdf. When the 3-parameter Weibull is used, the use of a lower limit is justified by the physical properties of the glass; instead, when the LTW is used, the lower limit is connected to the quality control to which glass panes are subjected at the end of the manufacturing process. In other words, the lower limit is determined by the maximum cracks size that can pass the quality control: since all plates with cracks deeper than a maximum allowable value are discarded, glass strength cannot drop below a minimum threshold value.Notwithstanding this recent modification, existing models are not yet satisfactory in case of glass plates produced with the float production process, patented in the ‘50s by Sir Alastair Pilkington [7, 8]. Briefly, in this industrial process, glass paste is poured on a bed of molten tin, forming a floating panel. Next, the glass plate is moved by steel rollers inside a lehr (i.e., a temperature-controlled kiln), where it is gradually cooled. The contact of the lower surface with the tin bath and the steel rollers (tin side) modifies the defectiveness scenario of the tin side that, at the end of the process, sensibly deviates from that of the surface exposed to the air (air side). More precisely, whereas the air side can be considered as pristine glass, the tin side can be assimilated to a predamaged surface. Because of this difference, the LTW can ensure an optimal fit (of the experimental data) only when the air side is under tension but not when it is the tin side to be under tensile stress [9].To close this gap in the literature, we propose a novel model based on the assumption that, during manufacturing, abrasion alters the original defect scenario of the pristine material, making unknown the distribution of the cracks’ length on the tin side. This assumption is supported by the experimental observation that glass strength, after abrasion, is no longer Weibull-shaped. The idea of a distortion that alters an original probability distribution is perfectly aligned with the change of variable theorem (CVT), sometime referred as “derivative distributions” too [10]. Indeed, given that (i) on the basis of the specialized literature, we expect to observe a Weibull pdf (for the pristine material) and (ii) the deformation affecting the data can be suitably represented through a function h(·), and the CVT provides the formula to generate the new pdf of the data, after the deformation.We will show that using the CVT-based protocol and without the need to make any assumptions on the distributions of the new cracks (which are known to be unobservable at least on an industrial level), experimental data are very well fitted, especially using a LTW and a 4-parameter logistic curve, as original pdf and transforming function, respectively. The selected transforming function is valuable also from a technical and engineering point of view, as its parameters have a clear physical interpretation, and can help practitioners in understanding the effect of abrasion on glass’ performance. Thus, our protocol provides a very suitable tool to statistically describe the strength of abraded glass, when abrasion is caused by the floating production process used to manufacture most of the glass panes applied in constructions.The work is organized as follows. Section2 provides a literature review about existing statistical models related to glass strength. Section 3 describes the sample at our disposal used to validate the fitting protocol. Section 4 deals with parameters estimation of the LTW distribution of the air side’s tensile strength, and it introduces the CVT. Section 5 contains our proposal of a function (the 4-parameter logistic) which suitably describes, in a statistical sense, the deformation caused by the abrasion process, together with the estimated pdf obtained through the CVT and its goodness of fit. Section 6 provides a physical and technical interpretation of the parameters of the model. The last two sections conclude and give directions for possible future researches.


## 2. Fitting Problem: A Brief Literature Review

As discussed in Introduction, the strong brittle behavior of glass panes is due to the unavoidable presence of flaws located on the external surfaces. These flaws are generally assumed to be “half-penny shaped” and they govern the mechanical response of the material. Collapse, in fact, occurs when the combination of the crack’s length and tensile opening stress reaches a critical condition, which rapidly leads to the macroscopic glass rupture. In this regard, we recall the existence of three different modes of crack’s propagation as follows:(1)
Opening mode (mode I) due to tensile stresses at right angle with crack axis(2)
Sliding mode (mode II) due to shear stresses parallel to the crack axis and perpendicular to the crack front(3)
Tearing mode (mode III) due to shear stresses parallel to both the crack axis and the crack frontIn case of glass panes, the contributions of both mode II and III are negligible, and collapse occurs when the combination of crack’s length and tensile opening stress reaches the critical condition given by(1)KIc=Υσπδ,where is the critical stress intensity factor (MPa m0.5) in opening mode I (KIc = 0.75 (MPa m0.5) for soda-lime silica glass [11]), Υ is the shape factor (Υ = 2.24/π for half-penny shaped crack), σ is the tensile stress (MPa) at right angle with the crack axis, and δ is the crack length (m) (the notation used in the remainder of the paper can be found in Appendix B).Equation (1) allows one to link the microscopic variable crack’s length with the macroscopic and easy to observe variable glass’ strength. As in [6, 12], if crack’s length is supposed to have a Paretian distribution with parameter α, glass’ strength can be modeled using a classical Weibull distribution. On this basis, from here on, we will focus on the random variable describing the strength, which is the main characteristic to be analyzed.Although the 2-parameter Weibull is certainly the most used model [13], it fails to properly fit the left-hand-side tail of the glass strength distribution [5, 14]. This fact represents the weak point of the model because, to be acceptable, structural components must ensure a very low failure probability. Hence, the left tail is the only part (of the glass strength distribution) to be of practical relevance in structural design.Starting from the consideration that glass elements must guarantee certain aesthetic and optical performances and that, consequently, glass plates must fulfill the factory quality control, a better model was proposed in [9]. Specifically, the authors assumed that, due to the rejection of panes with defects in transparency, cracks in the accepted panes cannot exceed a maximum allowable size. The effect of the upper truncation of the microcracks lengths introduces a lower bound in the glass resistance, whose distribution follows a LTW. In particular, the failure probability Pfσmax that a glass pane of area A breaks under a tensile stress lower or equal than σmax is given as(2)Pfσmax=1−exp−Akσmaxm−σ0mη0m,where m=2α−1 is the Weibull shape parameter and α is the Paretian parameter characterizing the population of microcrack lengths; η0=ΔAσk2α−1−σ02α−11/2α−1 is the Weibull scale parameter, with σ0 and σk corresponding to the stress obtained by plugging into equation (1), the cracks of maximal (δmax) and minimal size (δmin also known as physiological crack size), respectively. Both crack sizes can be determined following a graphical procedure, as detailed in [6]. Lastly, ∆A is a small portion of the area under stress. Specifically, ∆A is obtained by ideally dividing A into N representative elements, small enough to contain one microcrack at most, and k is a coefficient accounting for the deviation of the stress state from the equibiaxiality [15].As anticipated in Introduction, neither a left-truncated nor a generalized Weibull or any other classical statistical distribution provides a remarkable goodness of fit for the tin side of the glass panes [9]. Only recently, a significant improvement has been made in [6], where the problem of characterizing the abraded glass strength is faced by studying the effects of the superposition (mixture) of two different microcrack distributions. The first one characterizes the air-side defectiveness scenario, while the second one characterizes the cracks added by the abrasion process. Anyway, the goodness of fit provided by this statistic is barely acceptable. Indeed, performing a chi-squared goodness-of-fit test, the p value is around 6% and so, although the null hypothesis H0 (of proper fit) cannot be rejected, considering H0 as true, based on such a small p value, is rather heroic. The approach by [6] has another disadvantage, as it requires the preliminary choice of an appropriate pdf to model the length of the microcracks introduced by the abrasion process. The shape of such distribution is strongly affected by the events which cause the abrasion itself and so it is extremely hard to find a pdf that could be considered representative of almost all the abrasive processes. As we will explain in the following, the CVT makes it possible to bypass this problem, as no assumptions concerning the pdf of the unobservable microcrack length are needed. This is another reason that prompted us to exploit the potentialities offered by the CVT.


## 3. Data Description

In this study, we refer to the experimental campaign performed by the technical committee TC129/WG8 of the European Committee for Standardization (CEN) [16]. Glass plates, furnished by different European glass manufacturers, were tested under the coaxial double ring configuration with overpressure, according to the EN1288-2 standard [17]. A total of 740 valid failure stress values were obtained: 370 for the air-side strength and 370 for the tin side strength. This experimental campaign represents, albeit to our knowledge, the widest one performed on structural glass.As shown in Figure1, according to this test setup, specimens are supported by a ring and the load is transmitted through a second ring, coaxial with the supporting one.Figure 1
Basic diagram of the coaxial double ring test with overpressure (p) according to EN1288-2 [17].During the test, the overpressurep changes following the law specified by the standard [17], which correlates the overpressure p and the piston force in a way that should assure an equibiaxial state of stress in the core of the specimen (i.e., in the area bounded by the inner ring). In these conditions, the probability of finding the maximum tensile stress at right angle with the axis of the critical cracks is 100%. However, as shown in [18], the relationship among the overpressure p, the force applied by the tensiometer F, and the resulting tensile stress σTC129 given in the standard is not accurate. Furthermore, the real state of stress in the core of the specimen delimited by the inner ring is not equibiaxial. Thus, since all the relevant data are obtained from the value of the force F, but the correlation between the piston force and the tensile stress is not accurate, the stress values σTC129 recorded in the report [16] do not correspond to the maximum tensile stresses σmax acting within the specimens when the collapse occurs. For these reasons, rather than using the stress values σTC129 recorded in [16], we will refer to the corrected values reported in [18]. Specifically, in that work, via the finite element method (FEM) analysis, the correct correlation between F and σmax was found and, using the values of F associated with the stresses σTC129, the actual values of tensile stress σmax were obtained [18].


## 4. Air-Side Glass Strength and Change of Variable Theorem

### 4.1. Air-Side Distribution

At first, we verified that the LTW is a good choice to model the glass’ strength distribution on the air side. The corresponding density functiongXx is given by(3)gXx=A⋅k⋅mη0⋅xη0m−1⋅exp−Ak⋅xm−σ0mη0m.The corresponding cumulative distribution functionGX (x) is given by equation (2).In equation (3), A represents the area of the specimen bounded by the inner ring (in the testing configuration of Figure 1), and it is equal to 282600 mm (the same configuration will be used for all the other goodness-of-fit test mentioned in the following sections). k is the corrective coefficient of the effective area, and its value depends, mostly, on plate geometry, applied constraints, and type of loading. To be frank, k also depends on the value reached by the maximum tensile stress, which is produced by the actual combination of applied force F and overpressure p. Anyhow, if one neglects the nonlinear effects due to the membrane strains, the dependence of k on the amount of loading can be eliminated. By doing so, for the case at hand, a value of k = 0.45 was estimated via FEM analysis, as suggested in [9].All the other parameters ofgXx were graphically estimated on the air side, according to the method proposed in [6]. The following values were finally found:Shape factorm = 4.25Scale factorη0 = 1,462.85 (MPa mm2/α)Truncation valueσ0 = 40.2 [MPa]Using these values,gXx assures a very good fit of the experimental data, as graphically shown in Figure 2, where gXx is superimposed to the histogram of the glass strength on the air side.Figure 2
Observed histogram compared with the graphically estimated LTW (p value 78%).The quality of the obtained fit was ascertain performing aχ2-goodness-of-fit-test based on 20 equispaced intervals, that were grouped together on the tails of the distribution to avoid classes with less than 5 observations. In this way, a total of 15 classes were made (the same configuration will be used for all the other goodness-of-fit tests mentioned in the following sections). We also note that, since the LTW distribution has three parameters to be estimated (shape, scale, and truncation, whereas A and k are exogenous values), the test statistics is distributed, approximately, as χ15,3,12. The test returned a very high p value ≅ 78%, and so the null hypothesis H0 “the air-side strength sample is drawn from the LTW distribution provided bygXx” cannot be rejected.


### 4.2. Change of Variable Theorem

The crucial point is that one might expect to observe the same distribution of the glass strength both on the air and on the tin side. Conversely, as explained before, sample data on the tin side do not support the LTW. Since the LTW fits experimental data relative to the air-side strength, we claim the following:The original random variableX is deformed, due to the abrasion caused by the contact of the glass plate with the tin bath and the steel rollers (during the floating process)Data on the abraded side correspond to the distorted random variableY = h(X) with unknown pdf fYyConsequently, the main issue is no longer the choice of a suitable pdf to describeY, but rather the analytical description of the deforming function h(·). In light of this, the famous change of variable theorem [19, 20] perfectly mirrors the idea of analyzing the underlying deformation and describes “how areas are distorted by differentiable functions” [21]. In terms of pdf, the CVT assures that, given a transformed random variable Y = h(X) and provided that h(·) is continuous, differentiable, and monotonic, the pdf of Y can be obtained as follows:(4)fYy=ddyh−1y⋅gXh−1y.The original pdfgXx and the function h(·) are called, respectively, generating distribution and transforming function [22]. The choice of a proper generating function must be supported by solid technical arguments, widely accepted by the scientific community (as for the LTW in our case). If not so, it would not be possible to claim that data do not follow the expected generating function because of a distortion; lack of fit, in fact, could be simply related to a wrong choice of the original generating distribution.As far as we know, the CVT has never been used in the analysis of brittle materials, and, in our opinion, it represents a simple and powerful tool, which permits to statistically analyze data altered by abrasive processes.It is worth noting that abrasion is obviously not deterministic. Here, we are assuming thaty=hx+ϵ and that we can predict the average effect of abrasion on glasses with pristine strength x, whereas the randomness is described by the residuals ϵ. It is a reasonable assumption when the abrasion is due to the production process, as in our case. Here, the tin side strength is analyzed to show the potentialities of the approach, taking advantage of the large sample of data at our disposal. In the future, the same protocol will be suitable to be extended to interpret the variability of the strength of aged glass subject to several abrasive phenomena.


## 4.1. Air-Side Distribution

At first, we verified that the LTW is a good choice to model the glass’ strength distribution on the air side. The corresponding density functiongXx is given by(3)gXx=A⋅k⋅mη0⋅xη0m−1⋅exp−Ak⋅xm−σ0mη0m.The corresponding cumulative distribution functionGX (x) is given by equation (2).In equation (3), A represents the area of the specimen bounded by the inner ring (in the testing configuration of Figure 1), and it is equal to 282600 mm (the same configuration will be used for all the other goodness-of-fit test mentioned in the following sections). k is the corrective coefficient of the effective area, and its value depends, mostly, on plate geometry, applied constraints, and type of loading. To be frank, k also depends on the value reached by the maximum tensile stress, which is produced by the actual combination of applied force F and overpressure p. Anyhow, if one neglects the nonlinear effects due to the membrane strains, the dependence of k on the amount of loading can be eliminated. By doing so, for the case at hand, a value of k = 0.45 was estimated via FEM analysis, as suggested in [9].All the other parameters ofgXx were graphically estimated on the air side, according to the method proposed in [6]. The following values were finally found:Shape factorm = 4.25Scale factorη0 = 1,462.85 (MPa mm2/α)Truncation valueσ0 = 40.2 [MPa]Using these values,gXx assures a very good fit of the experimental data, as graphically shown in Figure 2, where gXx is superimposed to the histogram of the glass strength on the air side.Figure 2
Observed histogram compared with the graphically estimated LTW (p value 78%).The quality of the obtained fit was ascertain performing aχ2-goodness-of-fit-test based on 20 equispaced intervals, that were grouped together on the tails of the distribution to avoid classes with less than 5 observations. In this way, a total of 15 classes were made (the same configuration will be used for all the other goodness-of-fit tests mentioned in the following sections). We also note that, since the LTW distribution has three parameters to be estimated (shape, scale, and truncation, whereas A and k are exogenous values), the test statistics is distributed, approximately, as χ15,3,12. The test returned a very high p value ≅ 78%, and so the null hypothesis H0 “the air-side strength sample is drawn from the LTW distribution provided bygXx” cannot be rejected.


## 4.2. Change of Variable Theorem

The crucial point is that one might expect to observe the same distribution of the glass strength both on the air and on the tin side. Conversely, as explained before, sample data on the tin side do not support the LTW. Since the LTW fits experimental data relative to the air-side strength, we claim the following:The original random variableX is deformed, due to the abrasion caused by the contact of the glass plate with the tin bath and the steel rollers (during the floating process)Data on the abraded side correspond to the distorted random variableY = h(X) with unknown pdf fYyConsequently, the main issue is no longer the choice of a suitable pdf to describeY, but rather the analytical description of the deforming function h(·). In light of this, the famous change of variable theorem [19, 20] perfectly mirrors the idea of analyzing the underlying deformation and describes “how areas are distorted by differentiable functions” [21]. In terms of pdf, the CVT assures that, given a transformed random variable Y = h(X) and provided that h(·) is continuous, differentiable, and monotonic, the pdf of Y can be obtained as follows:(4)fYy=ddyh−1y⋅gXh−1y.The original pdfgXx and the function h(·) are called, respectively, generating distribution and transforming function [22]. The choice of a proper generating function must be supported by solid technical arguments, widely accepted by the scientific community (as for the LTW in our case). If not so, it would not be possible to claim that data do not follow the expected generating function because of a distortion; lack of fit, in fact, could be simply related to a wrong choice of the original generating distribution.As far as we know, the CVT has never been used in the analysis of brittle materials, and, in our opinion, it represents a simple and powerful tool, which permits to statistically analyze data altered by abrasive processes.It is worth noting that abrasion is obviously not deterministic. Here, we are assuming thaty=hx+ϵ and that we can predict the average effect of abrasion on glasses with pristine strength x, whereas the randomness is described by the residuals ϵ. It is a reasonable assumption when the abrasion is due to the production process, as in our case. Here, the tin side strength is analyzed to show the potentialities of the approach, taking advantage of the large sample of data at our disposal. In the future, the same protocol will be suitable to be extended to interpret the variability of the strength of aged glass subject to several abrasive phenomena.


## 5. Fitting the Tin-Side Data through the CVT

The next step consisted in the selection of the right form for the transforming function, aiming to efficiently describe the effect of abrasion on glass’ strength. We note that, for a proper selection of the functional form ofh(·), a trade-off should be reached between two different aspects of the problem. Indeed, if on the one hand the optimal choice should be physically motivated and supported by the knowledge of the abrasive phenomena, on the other one hand, the chosen function should have as few parameters as possible (a maximum of four is generally recommended [23]) to be easy to be estimated. We also note that, in case of simple transforming functions, such as polynomials, parameters can always be estimated through the least squares method, provided that the monotonicity constraint has been imposed.


### 5.1. Preliminary Considerations

To take into account the fact that theX values are distributed as a LTW, the regression analysis should be based on couples of quantiles xi,yii=1,…,n, as suggested in [24]. In detail, yi is the i-th of the n strength values observed on the tin side (370 in our case), sorted in the ascending order (i.e., the i-th empirical quantile). Similarly, xi is the corresponding theoretical quantile of the LTW function gXx, with parameters estimated as in Section 4.1 for the air side.This choice requires some relevant remarks. Firstly, to apply the CVT, we must suppose that the transforming function is monotonic. In our case, it seems fair to assume that glass plates which are more resistant than others will remain, on average, more resistant also after the abrasion. Hence, we are allowed to use quantiles to solve the regression model. Secondly, the use of quantiles makes it possible to analyze the way in which the abrasion acts on glasses located in different portions of the original Weibull distribution (for example, glasses corresponding to points on the left tail). Consequently, as it will be explained later, the use of quantiles gives us an insight about the order of magnitude of the cracks added by the same abrasion process.The resulting Q-Q plot is shown in Figure3, where theoretical quantiles of the original LTW (on the horizontal axis) are plotted versus the observed strength quantiles (on the vertical axis).Figure 3
Q-Q plot LTW theoretical quantiles versus the observed strength quantiles.For what was said in the previous sections, it is clear that a linear transformation would be useless. Indeed, ifh(x) = qx should fit well the couples of quantiles xi,yii=1,…,n, then, according to the CVT, the LTW should hold both on the air and the tin sides of the glass panes, as its functional form is preserved by multiplications by constant values. Nevertheless, the visual inspection of the Q-Q plot reveals that h(x) is almost linear for the highest quantiles. So, we decided to estimate the slope q of the regression line, limited to the 185 values higher than the median of the whole sample, equal to 62.33 MPa. The obtained regression line y = 0.712x is plotted as a solid line in Figure 3. For the sake of clarity, in the same figure, also the bisecting line y = x is shown as a dotted line. As it can be seen, the plot seems to confirm h(x) = qx as a suitable transformation for the subset of the higher values of the glass strength but not for the lower ones. A similar result was obtained in an alternative explorative analysis proposed in Appendix A (recall that glass strength is governed by microcracks on the surfaces, whose size distribution is a la Pareto, in case of unabraded surfaces; hence, the deformation of the strength distribution is a consequence of the deformation of the microcrack distribution; in Appendix A, the analysis of the deformation of the microcrack size distribution is recorded for a better comprehension of the phenomenon).


### 5.2. A Suitable Functional Form for h(·)

Based on the aforementioned exploratory results, in this section, we propose an ad hoc transformation, whose parameters have a relevant interpretation, in terms of the abrasion process. In particular, we suppose the transforming function has the following formula:(5)Y=hX=sX⋅X.Specifically, equation (5) makes explicit the ratio s(·) between the strength observed after (Y) and before (X) the abrasion. We expect that such ratio will be, on average, lower than 1, because glass’ strength should be reduced by abrasion. In addition, from the observation of Figure 3, we also formulate the following hypotheses:Abrasion during the production process adds some cracks on the glass pane, whose length is unobservable but lies in a given range of values which depends on the type of abrasive process.Very weak glasses (i.e., glass plates that were very weak also on the air side) cannot be further weakened through abrasion. In other words, their maximal crack is so large that the probability to add an even bigger crack is almost negligible. Hence, we pose thats(X) ˜ 1 for small values of X.Very strong glasses have very small cracks, which are negligible with respect to the newly added ones. Then, by reasonably assuming that new cracks are Pareto-like distributed, we can expect that the strength distribution after abrasion in the subset of highly strong glasses is againa la Weibull. Hence, we pose that, for large X, s(X) ⟶ q, where q is a constant value.Aiming to satisfy these requirements, we suggest using a sigmoid shape fors(·), that is a bounded, monotonic, and differentiable function, with two horizontal asymptotes. We also impose s(x) > 0 for x > 0, to make h(x) monotonic as required by the CVT.Now, to estimates(·) and h(·), some algebraic tricks are needed. Indeed, if s(·) is a sigmoid function (logistic, hyperbolic tangent, etc.), the inverse h−1(·) may not have an explicit formula. Nevertheless, if Y/X has a sigmoid trend (i.e., monotonic and bounded by two horizontal asymptotes), the reverse ratio X/Y must have the same behavior but with opposite monotonicity. Consequently, the previous hypotheses still hold, and we can suppose h−1Y=ξY⋅Y with ξ (Y) (called the transforming coefficient) having again a sigmoid shape. Finally, to satisfy the aforementioned requirements we choose, among all the possible sigmoid functions, the so-called 4-parameter logistic is given as(6)ξY=L4py=d+a−d1+y/cb,y>0,where a=L4p0, b is called Hill’s coefficient and it is related to the slope in the inflection point; consequently, it is a measure of the curve steepness, c is the halfway point between minimum and maximum, i.e., L4pc=0.5⋅a+d, and d is the height of the right horizontal asymptote.Also note thatL4py is monotonically increasing (respectively decreasing) if a < d (respectively a > d) and that d/dyL4pyy=0=0.To numerically estimate the coefficientsa, b, c, and d, we exploited the R’s package “dr4pl” [25, 26] which is able to fit logistic functions, using the observed values of Y, as explanatory variable, and the observed ratio X/Y, as dependent variable. The resulting estimated transforming function is(7)h−1y=ξy⋅y=1.41+1.08−1.411+y/50.4314.52⋅y.The outputs provided by theR package, together with the estimated confidence intervals (at a 95% confidence level), are listed in Table 1; Figure 4 displays the scatterplot of the ratio (ξy=x/y) versus the observed tin-side quantiles (y), together with the estimated logistic.Table 1
Summary of the estimated parameters of the transformed density functionfYy.
Left-truncated Weibull parameters+m (shape)η0 (scale)σ0 (truncation)4.251462.35 MPa mm1/α40.2 MPaL4p parameters++a (lower limit)b (slope)c (halfway point)d (upper limit)1.08 (0.92, 1.23)14.52 (2.49, 26.55)50.43 (46.85, 54.89)1.41 (1.39, 1.42)Additional parametersAk282600 mm20.45+Previously estimated through the graphical method. ++ Estimated with the R package dr4pl.Figure 4
Scatterplot of the ratio between the expected (X) and the observed (Y) strength.


### 5.3. Deformed pdf Fitting the Tin-Side Strength

PlugginggXx (equation (3)) and the just-now estimated h−1y (equation (7)) into equation (4), the resulting pdf fYy, namely, “Deformed Weibull”, is finally obtained. Table 1 contains a summary of all the parameters needed in the final formula of fYy. Since it is quite long to be explicated, for sake of simplicity, we propose the corresponding cdf FYy, which is more compact and can be immediately obtained by substituting x with h−1(y) in equation (2):(8)FYy=GXh−1y=1−exp−A⋅k⋅d+a−d/1+y/cb⋅y−σ0mη0m.Lastly, to evaluate the goodness-of-fit provided byfYy, we performed a χ2 test. The model has four parameters (because, in this case, A, k, m, η0, and σ0 are known), and so, under the null hypothesis H0 “the sample of observed strengths on the tin side is drawn fromfYy,” the test statistic is approximately distributed as a χ15,4,12. As in Section 4.1, to perform the test, the data were subdivided in 20 equispaced classes and grouped together to avoid near-to-zero frequencies, for a total of 15 classes. By doing so, a p value≅80% was found, and, consequently, H0 cannot be rejected. The quality of the fit is graphically shown in Figure 5, which indicates a very good adaptation also in the tails of the distribution.Figure 5
Histograms of the tin-side strength (MPa), compared with the deformed Weibull (p value 80%).It is evident that the CVT led to an exceptional goodness of fit with experimental data, when the tin side is the surface under tension. Conversely, nongeneralized Weibull nor other classical distribution is able to guarantee a good fit with the experimental tin-side strengths. We remind that the best model in literature is the bimodal truncated Weibull distribution, which furnishes ap value around 6% [5].


### 5.4. Robustness of the Model: Comparisons with Benchmark Distributions and Sensitivity Analysis

To further investigate the quality and the robustness of the proposed fitting protocol, we repeated the fitting procedure, both for the air and tin sides, using the benchmarks distributions cited in Section2 (i.e., 2- and 3-parameter Weibull and left-truncated Weibull).The obtained results are shown in Table2.Table 2
Fitting results of three benchmark distributions.
m (shape)η (scale)σ (location)p valueAir side2-parameter Weibull5.3667841.43-0.05383-parameter Weibull2.62974705.96390.7535Left-truncated Weibull4.25271462.8540.20.7832Tin side2-parameter Weibull6.5080415.27-1.24e-113-parameter Weibull2.60172847.31367.05e-09Left-truncated Weibull4.24741031.4835.712.06e-07Parameters estimated with the graphical approach;p value relative to a χ2 test with 15 classes.As it can be seen, none of the tested Weibull models is capable of fitting the distribution of the glass’ strength on the tin side; even the 2-parameter Weibull fails to properly fit the distribution on the air side. The same conclusions are graphically displayed in Figures6-7.Figure 6
Benchmark fitting distributions on the air side.Figure 7
Benchmark fitting distributions on the tin side.To conclude this section, we observe that the corrective coefficientk, which is a parameter of the LTW used as generating function (equation (3)), was estimated via FEM. Recall that its value depends on plate geometry, applied constraints, and type of loading. Furthermore, when geometric nonlinearities are of importance, it varies with the load level. However, under the configuration provided by standard EN 1288-2, geometrical nonlinearities are modest, i.e., k only slightly varies with the load level. In particular, for reliable values of failure stress values reached during the test, it is shown in ([5], Figure 6) that k could range in the following interval [0.4, 0.55]. So to ascertain the robustness of the obtained pdf, we checked the impact of k on the other parameters. Results are shown in Table 3.Table 3
Effects ofk on parameters of both LTW and 4p-logistic distributions.
Parameters of the LTWParameters of the 4p-logistickm/ηk/ησ/η1/η(m−1)Abcd0.4000.0029890.000281120.0274825.54E-111.08173714.5551750.412171.4077680.4250.0029460.000294460.0274825.29E-111.08171614.5555450.412311.4077660.4500.0029070.000307620.0274825.06E-111.0817414.5550950.412171.4077620.4750.0028700.000320610.0274824.86E-111.08172214.555150.411591.4077650.5000.0028360.000333440.0274824.67E-111.08173214.5554750.413021.4077640.5250.0028040.000346110.0274824.50E-111.0817414.5552250.412431.4077690.5500.0027730.000358650.0274824.34E-111.08167114.5552950.411141.407761Similar sensitivity analyses were not performed neither for the other parameters of the LTW nor for the 4-parameter logistic. All these parameters, in fact, were obtained through specific graphic or analytical procedures and are already optimized.


## 5.1. Preliminary Considerations

To take into account the fact that theX values are distributed as a LTW, the regression analysis should be based on couples of quantiles xi,yii=1,…,n, as suggested in [24]. In detail, yi is the i-th of the n strength values observed on the tin side (370 in our case), sorted in the ascending order (i.e., the i-th empirical quantile). Similarly, xi is the corresponding theoretical quantile of the LTW function gXx, with parameters estimated as in Section 4.1 for the air side.This choice requires some relevant remarks. Firstly, to apply the CVT, we must suppose that the transforming function is monotonic. In our case, it seems fair to assume that glass plates which are more resistant than others will remain, on average, more resistant also after the abrasion. Hence, we are allowed to use quantiles to solve the regression model. Secondly, the use of quantiles makes it possible to analyze the way in which the abrasion acts on glasses located in different portions of the original Weibull distribution (for example, glasses corresponding to points on the left tail). Consequently, as it will be explained later, the use of quantiles gives us an insight about the order of magnitude of the cracks added by the same abrasion process.The resulting Q-Q plot is shown in Figure3, where theoretical quantiles of the original LTW (on the horizontal axis) are plotted versus the observed strength quantiles (on the vertical axis).Figure 3
Q-Q plot LTW theoretical quantiles versus the observed strength quantiles.For what was said in the previous sections, it is clear that a linear transformation would be useless. Indeed, ifh(x) = qx should fit well the couples of quantiles xi,yii=1,…,n, then, according to the CVT, the LTW should hold both on the air and the tin sides of the glass panes, as its functional form is preserved by multiplications by constant values. Nevertheless, the visual inspection of the Q-Q plot reveals that h(x) is almost linear for the highest quantiles. So, we decided to estimate the slope q of the regression line, limited to the 185 values higher than the median of the whole sample, equal to 62.33 MPa. The obtained regression line y = 0.712x is plotted as a solid line in Figure 3. For the sake of clarity, in the same figure, also the bisecting line y = x is shown as a dotted line. As it can be seen, the plot seems to confirm h(x) = qx as a suitable transformation for the subset of the higher values of the glass strength but not for the lower ones. A similar result was obtained in an alternative explorative analysis proposed in Appendix A (recall that glass strength is governed by microcracks on the surfaces, whose size distribution is a la Pareto, in case of unabraded surfaces; hence, the deformation of the strength distribution is a consequence of the deformation of the microcrack distribution; in Appendix A, the analysis of the deformation of the microcrack size distribution is recorded for a better comprehension of the phenomenon).


## 5.2. A Suitable Functional Form for h(·)

Based on the aforementioned exploratory results, in this section, we propose an ad hoc transformation, whose parameters have a relevant interpretation, in terms of the abrasion process. In particular, we suppose the transforming function has the following formula:(5)Y=hX=sX⋅X.Specifically, equation (5) makes explicit the ratio s(·) between the strength observed after (Y) and before (X) the abrasion. We expect that such ratio will be, on average, lower than 1, because glass’ strength should be reduced by abrasion. In addition, from the observation of Figure 3, we also formulate the following hypotheses:Abrasion during the production process adds some cracks on the glass pane, whose length is unobservable but lies in a given range of values which depends on the type of abrasive process.Very weak glasses (i.e., glass plates that were very weak also on the air side) cannot be further weakened through abrasion. In other words, their maximal crack is so large that the probability to add an even bigger crack is almost negligible. Hence, we pose thats(X) ˜ 1 for small values of X.Very strong glasses have very small cracks, which are negligible with respect to the newly added ones. Then, by reasonably assuming that new cracks are Pareto-like distributed, we can expect that the strength distribution after abrasion in the subset of highly strong glasses is againa la Weibull. Hence, we pose that, for large X, s(X) ⟶ q, where q is a constant value.Aiming to satisfy these requirements, we suggest using a sigmoid shape fors(·), that is a bounded, monotonic, and differentiable function, with two horizontal asymptotes. We also impose s(x) > 0 for x > 0, to make h(x) monotonic as required by the CVT.Now, to estimates(·) and h(·), some algebraic tricks are needed. Indeed, if s(·) is a sigmoid function (logistic, hyperbolic tangent, etc.), the inverse h−1(·) may not have an explicit formula. Nevertheless, if Y/X has a sigmoid trend (i.e., monotonic and bounded by two horizontal asymptotes), the reverse ratio X/Y must have the same behavior but with opposite monotonicity. Consequently, the previous hypotheses still hold, and we can suppose h−1Y=ξY⋅Y with ξ (Y) (called the transforming coefficient) having again a sigmoid shape. Finally, to satisfy the aforementioned requirements we choose, among all the possible sigmoid functions, the so-called 4-parameter logistic is given as(6)ξY=L4py=d+a−d1+y/cb,y>0,where a=L4p0, b is called Hill’s coefficient and it is related to the slope in the inflection point; consequently, it is a measure of the curve steepness, c is the halfway point between minimum and maximum, i.e., L4pc=0.5⋅a+d, and d is the height of the right horizontal asymptote.Also note thatL4py is monotonically increasing (respectively decreasing) if a < d (respectively a > d) and that d/dyL4pyy=0=0.To numerically estimate the coefficientsa, b, c, and d, we exploited the R’s package “dr4pl” [25, 26] which is able to fit logistic functions, using the observed values of Y, as explanatory variable, and the observed ratio X/Y, as dependent variable. The resulting estimated transforming function is(7)h−1y=ξy⋅y=1.41+1.08−1.411+y/50.4314.52⋅y.The outputs provided by theR package, together with the estimated confidence intervals (at a 95% confidence level), are listed in Table 1; Figure 4 displays the scatterplot of the ratio (ξy=x/y) versus the observed tin-side quantiles (y), together with the estimated logistic.Table 1
Summary of the estimated parameters of the transformed density functionfYy.
Left-truncated Weibull parameters+m (shape)η0 (scale)σ0 (truncation)4.251462.35 MPa mm1/α40.2 MPaL4p parameters++a (lower limit)b (slope)c (halfway point)d (upper limit)1.08 (0.92, 1.23)14.52 (2.49, 26.55)50.43 (46.85, 54.89)1.41 (1.39, 1.42)Additional parametersAk282600 mm20.45+Previously estimated through the graphical method. ++ Estimated with the R package dr4pl.Figure 4
Scatterplot of the ratio between the expected (X) and the observed (Y) strength.


## 5.3. Deformed pdf Fitting the Tin-Side Strength

PlugginggXx (equation (3)) and the just-now estimated h−1y (equation (7)) into equation (4), the resulting pdf fYy, namely, “Deformed Weibull”, is finally obtained. Table 1 contains a summary of all the parameters needed in the final formula of fYy. Since it is quite long to be explicated, for sake of simplicity, we propose the corresponding cdf FYy, which is more compact and can be immediately obtained by substituting x with h−1(y) in equation (2):(8)FYy=GXh−1y=1−exp−A⋅k⋅d+a−d/1+y/cb⋅y−σ0mη0m.Lastly, to evaluate the goodness-of-fit provided byfYy, we performed a χ2 test. The model has four parameters (because, in this case, A, k, m, η0, and σ0 are known), and so, under the null hypothesis H0 “the sample of observed strengths on the tin side is drawn fromfYy,” the test statistic is approximately distributed as a χ15,4,12. As in Section 4.1, to perform the test, the data were subdivided in 20 equispaced classes and grouped together to avoid near-to-zero frequencies, for a total of 15 classes. By doing so, a p value≅80% was found, and, consequently, H0 cannot be rejected. The quality of the fit is graphically shown in Figure 5, which indicates a very good adaptation also in the tails of the distribution.Figure 5
Histograms of the tin-side strength (MPa), compared with the deformed Weibull (p value 80%).It is evident that the CVT led to an exceptional goodness of fit with experimental data, when the tin side is the surface under tension. Conversely, nongeneralized Weibull nor other classical distribution is able to guarantee a good fit with the experimental tin-side strengths. We remind that the best model in literature is the bimodal truncated Weibull distribution, which furnishes ap value around 6% [5].


## 5.4. Robustness of the Model: Comparisons with Benchmark Distributions and Sensitivity Analysis

To further investigate the quality and the robustness of the proposed fitting protocol, we repeated the fitting procedure, both for the air and tin sides, using the benchmarks distributions cited in Section2 (i.e., 2- and 3-parameter Weibull and left-truncated Weibull).The obtained results are shown in Table2.Table 2
Fitting results of three benchmark distributions.
m (shape)η (scale)σ (location)p valueAir side2-parameter Weibull5.3667841.43-0.05383-parameter Weibull2.62974705.96390.7535Left-truncated Weibull4.25271462.8540.20.7832Tin side2-parameter Weibull6.5080415.27-1.24e-113-parameter Weibull2.60172847.31367.05e-09Left-truncated Weibull4.24741031.4835.712.06e-07Parameters estimated with the graphical approach;p value relative to a χ2 test with 15 classes.As it can be seen, none of the tested Weibull models is capable of fitting the distribution of the glass’ strength on the tin side; even the 2-parameter Weibull fails to properly fit the distribution on the air side. The same conclusions are graphically displayed in Figures6-7.Figure 6
Benchmark fitting distributions on the air side.Figure 7
Benchmark fitting distributions on the tin side.To conclude this section, we observe that the corrective coefficientk, which is a parameter of the LTW used as generating function (equation (3)), was estimated via FEM. Recall that its value depends on plate geometry, applied constraints, and type of loading. Furthermore, when geometric nonlinearities are of importance, it varies with the load level. However, under the configuration provided by standard EN 1288-2, geometrical nonlinearities are modest, i.e., k only slightly varies with the load level. In particular, for reliable values of failure stress values reached during the test, it is shown in ([5], Figure 6) that k could range in the following interval [0.4, 0.55]. So to ascertain the robustness of the obtained pdf, we checked the impact of k on the other parameters. Results are shown in Table 3.Table 3
Effects ofk on parameters of both LTW and 4p-logistic distributions.
Parameters of the LTWParameters of the 4p-logistickm/ηk/ησ/η1/η(m−1)Abcd0.4000.0029890.000281120.0274825.54E-111.08173714.5551750.412171.4077680.4250.0029460.000294460.0274825.29E-111.08171614.5555450.412311.4077660.4500.0029070.000307620.0274825.06E-111.0817414.5550950.412171.4077620.4750.0028700.000320610.0274824.86E-111.08172214.555150.411591.4077650.5000.0028360.000333440.0274824.67E-111.08173214.5554750.413021.4077640.5250.0028040.000346110.0274824.50E-111.0817414.5552250.412431.4077690.5500.0027730.000358650.0274824.34E-111.08167114.5552950.411141.407761Similar sensitivity analyses were not performed neither for the other parameters of the LTW nor for the 4-parameter logistic. All these parameters, in fact, were obtained through specific graphic or analytical procedures and are already optimized.


## 6. Role of the 4p-Logistic Parameters

A very relevant fact is that each one of the four parameters of the L4p function plays a specific role in describing the abrasion phenomenon, which caused the deformation. As explained before, due to practical needs, we estimated the parameters of the inverse relationship X=h−1Y rather than Y=hX. Nevertheless, the physical meaning of each parameter can be easily reconstructed, as detailed below.First, we observe thata = 1.08 is statistically equivalent to 1, as confirmed by the 95% confidence interval (0.92, 1.23) provided in Table 1. This confirms that, at the lowest quantiles, where L4py≈a, the glass’ strength is almost unaffected by the transformation. The ratio Y/X of the strength after and before the abrasion is, in fact, very close to 1 (Y/X 1/a = 0.92), and so the additional weakening of glass panes there were already “intrinsically weak” is very small. From a physical perspective, this means that the contact with the tin bath and the steel rollers does not generate the very large cracks that are associated with the observed minimum values of the strength.On the other one side, the parameterd corresponds to the right horizontal asymptote for the ratio X/Y. Then, limited to the subset of the “very strong” glasses, we claim that the ratio Y/X 1/d = 0.712, which coincides with the value q, the slope of the linear regression shown in Figure 3. We recall that the linear regression was limited to the 185 specimens with strength on the tin side higher than the median value of 62 MPa, a fact that corresponds to the visual analysis of Figure 4, which shows that the S-shaped curve almost coincides with its horizontal asymptote exactly after this value of the strength. Also, this fact indicates that most of the cracks generated by abrasion are of modest size, but they are larger, on average, than those typical of the strongest specimens. The more severe is the abrasion process, and the highest is the expected value for parameter d.Parameterc is the point where the logistic function reaches the halfway value (a + d)/2 and parameter b is the steepness of the curve. It is of interest to observe that the logistic function reaches the one-fourth-way point (a + d)/4 at abscissa c1/3b and the three-fourth-way point 3 (a + d)/4 at abscissa c3b. The distance between these two abscissas decreases when b increases or c decreases. Thus, b and c indicate how quickly the logistic function passes from its minimum to its maximum and vice versa. Together, these two parameters give a measure of the dispersion of the cracks generated by the abrasion process, as they somehow define the intermediate range of glasses (neither very weak nor very strong) whose pre-existing cracks are comparable with the newly ones added through abrasion. Thus, b and c furnish information about the magnitude of the unobservable cracks caused by the abrasion. One can expect that, if abrasion generates homogeneous cracks (same size), the S-shaped function will approximate a step function. In that case, indeed, the transition from the distribution of strengths governed by the pre-existing cracks to the one governed by the newly added ones is sudden, and so b is expected to be high. On the other hand, if the abrasion process is not homogeneous and new cracks lengths cover a wide range of values overlapping the pre-existing ones, b is expected to be small and hence the transition from the original left-truncated Weibull to the distorted one is less neat and sharp.


## 7. Additional Technical Considerations

As noted in Introduction, to properly account for the effect of aging and wear on load-bearing structures, engineers make frequent use of corrective factors that can be found in international standards and/or in technical guidelines. Unfortunately, given the very recent introduction of structural glass components, practical knowledge concerning the behavior of these structures is still limited and the available corrective factors cannot be considered very reliable. Consequently, a statistical model capable to describe the variability of aged glass’ strength would be of paramount importance, as it would allow to statistically recalibrate the abovementioned corrective factors. In this regard, we noted that many scientific works have shown that the tensile strength of glass panes is distributed as a Weibull distribution (specifically as a left-truncated Weibull). However, this holds only for pristine material and not for aged glass. So we have proposed a novel protocol based on the use of the notorious change of variable theorem, that allows to determine the way in which the original distribution function of pristine glass is deformed due to the effect of abrasion. To test our protocol, we made use of a large sample of glass plates produced with the commonly used “float-manufacturing” process. The rational among this choice can be traced in the fact that, during manufacturing, the surface in contact with the tin bath and the steel rollers has actually undergone an abrasion process. Obtained results are encouraging, as the fitting of the tin-side strengths is excellent and far superior than that obtainable with other statistical models proposed in the literature. Hence, we expect that the protocol could be applied with similar success also to describe other abrasion processes, to which structural glass could be subject during its lifetime (such as the effect of sandstorms on buildings located in desertic regions). Certainly, to prove this conjecture, a properly designed experimental campaign is desirable, so as to collect data on different specimens subject to alternative abrasion processes and at different severity levels. This could be a topic for future research activities.


## 8. Conclusions

In this paper, we introduced a novel fitting protocol (based on the change of variable theorem) that makes it possible to statistically describe the tensile strength of aged and/or abraded glass. Our protocol is based on the assumption that the abrasion phenomenon acts as a deforming process on the statistical distribution of the depth of the existing surface cracks and, consequently, on the distribution of the glass’ strength. From technical literature, we knew that glass’ resistance is distributed as a left-truncated Weibull (LTW) and so, our main contribution was the definition of a continuous, differentiable, and monotonic function which suitably describes the way in which the original LTW is deformed through abrasion. Specifically, we proposed a 4-parameter logistic function, whose sigmoidal shape led to a very good fit of the tin-side strength distribution of a very large sample of glass plates used as benchmark.The use of the 4-parameter logistic function has an additional merit, as its parameters have an immediate physical meaning. Indeed, from the analysis of the obtained coefficients, we found that the population of abraded glass’ strength preserves a LTW shape for the higher quantiles, whereas it deviates from this function for the lower ones. A reasonable interpretation of this finding is that the stronger elements, characterized by small pre-existing cracks, are significantly weakened through abrasion, since the original cracks are negligible with respect to the newly added ones. Hence, new cracks generated by the abrasion process govern the higher quantiles of the strength distribution. On the other one hand, it is very unlikely that cracks generated by the abrasion process could be bigger or at least comparable with the largest one, already present in the weakest specimens.



---

*Source: 8519426-2020-07-29.xml*