An analytical solution for pre-crack behaviour of laminated glass under blast loading

Laminated glazing is often employed to minimise damage and injuries during blast events. In this work, the von Karman theory for large deﬂections of plates was used to simulate the effect of large explosions on laminated glazing. Linear material properties were assumed for both the glass and Polyvinyl Butyral layers. The glass and PVB layers were assumed to act fully compositely during the pre-crack phase of the deformation. A higher order deﬂection function was employed to represent the complex deformed shape observed in DIC blast test data collected by Hooper et al. (2012). The deﬂection results showed that the method developed could produce accurate estimates of the glazing deformation history during a blast event. The analytical solution was also used to compute the reaction forces acting on the window supports, which were found to be of a similar magnitude to those calculated from experimental data. In addition, crack densities were predicted, which were found to follow a pattern similar to those seen in blast experiments. The analytical approach developed is valuable for risk assessment engineers and façade designers who much prefer analytically based models over full-scale FE analysis, as FEA is often too time consuming for design assessments. 2016 The


Introduction
The façade system of buildings is of high importance for designers attempting to improve the security of structures against accidental or malicious explosions. The outer layer of the affected buildings needs to prevent blast waves from penetrating the building interior, protecting the occupants and internal fittings. Windows often represent particularly weak areas in this context. Traditional annealed glass is not well suited to protect from blast shock waves, as it tends to fracture early, projecting fragments inside and outside the building and allowing further blast pressure waves to penetrate the building envelope. Laminated glass has been shown to significantly increase the protection of structures [2]. This composite material is formed by two or more layers of glass interposed with layers of Polyvinyl Butyral (PVB), a polymer membrane with good optical properties. When loaded with the pressure wave, the glass layers will tend to fracture. Following this, the PVB membrane is able to stretch, absorbing a significant quantity of incident energy. Additionally, the glass fragments will generally remain attached to the interlayer, preventing them from being flown into the building space.
However, this composite material is difficult to account for in design, as its behaviour is complex. Additionally, the impulsive nature of the loading and the large deflections caused invalidate several simple analytical techniques. To obviate this, single degree of freedom approximations are often used, employing several empirical constants to obtain the required design data [2,3]. To obtain such constants and to study the effect of such loadings on window panes, blast experimental programmes have been carried out in the past [4]. Hooper et al. [1] performed tests where Digital Image Correlation (DIC) techniques were also used to collect deformation data throughout the extent of the glazing pane. These data was used subsequently by Del Linz et al. [5] to calculate the reaction forces acting on the supports and therefore provide further data to determine empirical constants to be used in design. Several other authors [6][7][8][9][10] also developed finite element analysis (FEA) models to perform analyses of the system. Whilst the results obtained are often realistic, the models can become very complex, requiring significant expertise and computing power to be utilised.
Wei and Dharani [11] developed an analytical model to predict the large deformations of the glazing before the failure of the glass layers, assuming a first order cosine wave to represent the deflected shape. They used this model to reproduce results from several experiments, which generally utilised relatively small charges. The method was then also extended to calculate likely crack densities in different areas of the pane [12]. Whilst their results agreed well with the experimental data for low explosive quantities, the solution proposed was not able to account for the behaviour shown by the DIC results for heavier blast loadings. In these, the measured deflected shape of the window indicated that a single term solution would not be able to represent the deflections accurately, as higher order cause the deflected shape to differ from a simple sinusoidal deflection in the early phase of the loading [1]. Therefore, in this paper, the analytical solution was expanded to account for these effects. Both the deflection and the Airy's stress function were considered as summations of several terms, and a general solution for any number of terms was derived. This was employed to produce the required differential equations, which were solved numerically. The deflected shapes were then compared with the experimental results. Additionally, the results were used to calculate the stresses acting on the window pane, and from these reaction forces and likely crack densities were found. These again were compared with experimental results available, with the aim of validating the proposed method.

Method
Three of the blast tests (Tests 1-3) performed by Hooper on laminated glass were considered for this analysis [1]. One of these tests (Test 1) used 12.8 kg of C4 explosive, (15 kg TNT equivalent), whilst two more tests (Test 2 and Test 3) employed 25.6 kg of C4 explosive, (30 kg TNT equivalent). The charge weights and explosive stand-offs are summarised in Table 1. All the glazing samples were 1.5 Â 1.2 m in planar dimensions and were composed of two 3 mm thick layers of annealed glass and a 1.52 mm PVB layer. Data was collected using pressure gauges, strain gauges and 3D DIC. This last technique allowed the measurements of full field 3D deflections throughout the panes. The experimental set up is shown in Fig. 1, whilst typical DIC results are shown in Fig. 2.
These data collected by Hooper et al. [1] were employed to validate an analytical solution of the plate deformation developed using the Von Karaman theory for large deformations, a commonly used method to represent these phenomena [11,13,14]. The plate was assumed to be thin, therefore through thickness effects were ignored in the analysis. Using this theory, however, required some assumptions both regarding the material properties of the panel and the loading. These will be considered first.

Material properties
Laminated glass is composed of two or more layers of glass interposed with layers of Polyvinyl Butyral, a rubbery polymer. Whilst any type of glass could be used, in this analysis annealed glass was assumed, as this was the variety used in the blast experiments.
Different levels of composite deformation theory or action between the composite components were also presumed. How-ever, as in previous studies full composite action was considered a reasonable assumption previous to the glass failure point [1,12]. Therefore, the rule of mixtures could be used to estimate the stiffness and Poisson's ratio of the system. The formulae used were: where E g and E p were the Young's modulus of the glass and PVB respectively, m g and m p were the Poisson's ratios and h g and h p were the thicknesses of each layer. The mass per unit area of the section could be found through: where q g and q p were the density of glass and PVB respectively.
Therefore, to apply the formulae above, it was necessary to assume appropriate material properties for each material. Whilst the glass layers would be normally considered as acting in a linear elastic manner previous to failure, the PVB membrane material was significantly more complex to model. Authors performed several experimental studies on this material, which showed significant non linearity and rate dependency in its behaviour [15,16]. Including a material model to accurately represent this complexity in the analytical structural model would have proven impractical. However, before the glass failure, the strains of the PVB membrane were small and the stresses significantly lesser than those in the glass layers, due to the much higher stiffness of the outer material. Because of this, previous authors normally considered the PVB as linear elastic for the purpose of pre-crack analysis [5]. The same assumption was made here, and the relevant material properties were obtained from the work of Hooper et al. [1]. The material properties assumed for the analysis are summarised in Table 2.
Part of the analysis also covered the analysis of the glass fragment density in the windows area and the calculation of reaction forces at the edges. To decide at which time point the crack density was to be calculated and the deflection and reaction calculations stopped, a failure condition for the glass plies had to be assumed. Glass failure is difficult to predict accurately, as it strongly depends on the presence of local random imperfections. Previous authors have used tensile stress levels of 80-100 MPa as the cracking limit for glass subject to explosions [1,12]. In this case an overall limit of 100 MPa was assumed. The resulting failure times were then compared with the experimental results, where the failure of the window could be determine within an error of ±1 ms from the DIC data. Once the failure time was decided it was used as described above, as it was also supposed that, as seen in the experiments, most of the crack surfaces would have been formed in less than 1 ms from the initial crack formation.

Blast properties
The pressure time history of blast pressure waves is generally composed of an initial near instantaneous rise, followed by an approximately exponential decay which eventually lowers the pressure to below atmospheric. Henceforth, a negative phase takes place, where the pressure on obstacles is reversed. In this analysis,  the blast pressures were included using the blast curve utilised in several previous studies [11,13,14,17]. The equation representing the blast wave was [18]: where p 0 is the initial pressure peak above the atmospheric pressure (overpressure), t is the time, t d is the time duration of the positive phase and a is a decay constant.
As the reflected pressure information was not directly available for all the blast tests considered, for this analysis the computational fluid dynamics results also used by Hooper et al. [1] where used. Specifically, the peak pressure p 0 and impulse data were used to fit the parameters of the blast wave equation for each case. The final constants used are shown in Table 3.

Analytical method
The analysis method used needed to allow for the large deflections which took place during the loading. This meant that membrane actions were likely to be present in the deformed system even before the glass failure and could not be ignored in the solution. The von Karman large deflection theory could take this into consideration and therefore was used in this study. The basic equations describing the motion were [11]: Dr 4 w þ Mw ;tt À hðw ;xx F ;yy þ w ;yy F ;xx À w ;xy F ;xy Þ where h was the total thickness of the sample, w was the out of plane deflection, F was the Airy's stress function and M was the mass of the section. In this solution, it was assumed that the glass would be simply supported along all four edges, without allowing any movement in the out-of-plane or in-plane direction. Therefore, assuming that the coordinates 0,0 were at the centre of the panel, a was the panel length and b was the panel width, the conditions reduced to: Additionally, an initial condition of zero deflections and rotations throughout the panel at t = 0 was utilised in the solution.
Dharani and other authors [11,13,14,17] assumed a deflected shape given by a combination of cosine waves. However, most of the authors, except for Türkmen and Mecitoglu, assumed a single term for their series, arguing this would be sufficient to represent the pre-cracking deformations of the pane. Whilst this assumption was appropriate for lower intensity blasts, the DIC results showed that such a deflection shape was not suitable for more intense loading situations. In these cases the curvature was not uniformly distributed throughout the panel, being instead concentrated in bands along the edges, as can be seen in Fig. 2. This was especially the case in the earlier stages of the deformation. To model this appropriately, it was decided to utilise more terms in the series. The equation was therefore assumed to be a summation: where u mn (t) were time dependent functions to be determined and m and n were odd positive integers. Eq. (8) could then be substituted in Eq. (5) giving: Therefore a summation for F of a form similar to the function used by previous authors [11] could be assumed: where f mn and k mn were coefficients to be determined. Differentiating Eq. (10) and substituting into Eq. (9) produced: This could be manipulated to give: Eq. (6) could be changed into: Lðw; FÞ ¼ Dr 4 w þ Mw ;tt À hðw ;xx F ;yy þ w ;yy F ;xx À w ;xy F ;xy Þ The Galerkin method could then be applied to Eq. (15): Eqs. (14) and (8) could be substituted in Eqs. (15) and (16). When proceeding with the derivation of the differential equation, it was possible to use the orthogonality of the cosine terms in the following situation: which, for odd positive n and m, was equal to 0 unless n = m. Using this result and following some simplifications, the following system of differential equations was obtained: Table 3 The blast curve constants used in the analytical calculations and the blast impulse for each case considered. where d represents the Dirac function. Eq. (18) was solved using the numerical differential equation solver capabilities of Matlab [19]. It was decided to use m and n up to 7 as this combination produced results which did not change significantly when further increasing the number of terms.

Reactions
The out of plane reactions imposed on the window edge are of prime importance to the design of the window system. Whilst most of the blast energy is normally assumed to be absorbed in the post-crack stage of the deformation, previous results suggested that significant reactions could be imposed in the pre-crack stage [5]. The experimental data used to obtain the reaction forces though presented significant noise, preventing an accurate calculation of the forces involved. Therefore it was considered useful to calculate these directly from the results of this analytical analysis.
It was assumed that both the bending and the membrane mechanisms produced significant forces in the pre-crack phase and needed to be considered in the analysis. A diagram showing the relevant forces is shown in Fig. 3.
The bending stresses were given by: where r bx , r by and s bxy are respectively the stresses in the x and y direction and planar shear. z is the position in the plate thickness and was assumed to be 0 at the centre. Eqs. (19)-(21) represent an approximation, as the stresses will not be constant at the interface between the glass and the PVB, however the level of errors introduced was considered acceptable for this analysis. The membrane stresses are instead given by: where r mx , r my were the membrane stresses in the x and y directions.
Finding the contribution to the reactions of the membrane forces required firstly calculating the angle of the glass at the edges, as it was assumed that these reactions would be acting parallel to the glass surface. This angle was calculated at all locations using the deflection data. The results were then used to calculate the out-of-plane component of the membrane stress. The stresses at each point were then integrated to find the reaction force. An area given by the glass thickness and the calculation points spacing was used for each datum.
Reactions due to the bending mechanism were taken to be equal to the shear force acting at the edge locations. In turns, this shear force was assumed to be equal to the differential of the bending moment. Therefore, to calculate this component of the reactions, firstly the bending moment was found from the bending stresses assuming a symmetrical distribution across the sample thickness. This bending moment was then differentiated in the x and y directions to find the shear forces at all points. The values at the edges were considered to be the out of plane reactions due to bending.
The membrane and bending reactions found in this manner where then integrated and compared with the experimental values.

Crack spacing
In all the cases considered the glass fractured during the loading. As the location of these cracks could help to assess the behaviour and residual capacity of the system after the failure of the glass, the stress results available were used to calculate the likely crack density on the window area. This was compared with descriptions of the experimental outcomes.
To achieve this, a method inspired by that used by Wei and Dharani [12] was applied. The energy balance of the glazing system was given by: where U is the total energy, U 0 is the strain energy, U a is the drop in energy due to a crack being formed and U c is the surface energy required to create a crack. When a crack is forming Eq. (24) reduced to: The surface energy required to form a crack was given by: where l is the length of crack, h g is the glass thickness and c s is 3.9 J/m 2 as given in [20]. The strain energy available in a facet was: ij / rs m 2 i 2 dðs À n À 2jÞ þdðs þ n À 2jÞ þ dðs À n þ 2jÞ dðr À mÞ where Dx and Dy are the dimension of each facet. If these two energies are equated, the value of l can be found. This represented the average length of cracks present in one of the calculation facets and was used as a measure of the cracking density. As the crack density was not measured during the blast tests, the results of this analysis could only be compared qualitatively with the description and images of typical results.

Results
The formulae derived above were used to calculate the deflection of the window pane. A quarter of the area was used for the calculations as, due to symmetry, the results would have been the same in all other quadrants. The calculated deflections were compared with the experimental DIC results to ensure that both the shape and magnitude of the reactions were comparable. Typical results along a cut located as shown in Fig. 4 are shown in Figs. 5, 6 and 7. The results for an equivalent cut taken in the opposite direction were similar. The analytical results managed to predict the general behaviour of the glass, generally showing a similar deflected shape as seen in the experimental results.
Some differences could be seen at the edges of the window panes. These were due to the idealisation of the boundary conditions. In reality the supporting steel members showed some deflections, though this was not measured directly. It can though be seen in the experimental results that the deflections do not always start from 0 at the glazing edges. As this was not included in the analytical model, the edge locations and central deflections are often smaller than found in the experiments. These differences are though limited to approximately 15% in the worst locations away from the edges. Additionally, DIC data was easily lost near the supports and part of the experimental data had to be interpolated, as described by Del Linz et al. [5]. This would have also introduced some irregularities.
The results were then used to calculate the stresses across the window pane. As quarter symmetry was used, the results plots show the data for an area as shown in Fig. 8. Fig. 9 shows typical results for Test 2 at 2.25 ms. The results at this time peaked at approximately 100 MPa, suggesting the window was likely to fail at approximately this moment. Whilst this was approximately 0.5 ms earlier than assumed from the experiments, such a difference is compatible with the uncertainty in the assessment of the experimental results. The failure time for the other blasts varied, being affected by the explosion energy. In Test 1 the failure time was higher, approximately 2.5 ms, whilst in test 3 due to higher blast intensity the failure took place at 1.75 ms.
Following this, the reaction forces were considered. The bending moments across the slab were firstly calculated. Typical results are shown in Fig. 10. As could be expected from the deflected shape and stress calculations, the peak moments were located in a band approximately 1/3 of the distance from the supports to the centre of the pane. This was also the area with the highest curvature in the deflections.
This data was differentiated and the shear forces were found along each corner. These were then added with the membrane forces, providing an overall estimate of the out-of-plane reactions. The results for the three tests are shown in Fig. 11. The analytical  results showed a similar magnitude as those seen in the experimental results, though with some deviations. Test 2 showed higher reactions in the analytical case, whilst in Tests 1 and 3 the analytical and experimental data showed a better agreement. The differences though could be due to the paucity of experimental time steps which could be used, only three except for test 1.
The stress data was also used to calculate likely crack densities in the panes. Typical results are shown in Fig. 12, together with an image from a typical experiment. The data showed a much higher density of cracks near the corner, followed by bands approximately 1/3 of the way between the edges and the centre of the window, with very few cracks present in the centre of the pane. This was similar to the previously observed stress concentrations, and was comparable to the distribution suggested by experimental observations [21].

Discussion
The results of the analysis showed that considering more than one term in the deflection function was essential to capture the pane behaviour in the pre-crack phase for the large blasts considered. Whilst lower blast intensities are likely to excite simpler deflection modes which could be adequately represented by the first term in the cosine series, the more extreme excitation of the structure in these heavier loading cases caused the glass to deform in a non-uniform manner. A consequence was that the highest bending moments were generated away from the glass centre. This affected the stress distribution throughout the pane, causing significant differences both in the reaction forces and cracking pattern, both aspects of direct interest to glass façade designers.
The calculated deflected shapes agreed with the experimental records. Whilst the absolute quality varied between the tests, the general behaviour of the pane was captured in all cases, with some showing very accurate predictions. Improvements could be made by refining the boundary conditions modelling. It is likely that in reality the edges would show a small degree of rigidity, which is not represented here. The results showed that a fully pinned model represented the results adequately, however introducing some moment resistance might improve the predictions. More   9. Stresses at the top surface across the window pane for Test 2 at 2.25 ms. The area covered by the plot is shown in Fig. 8.   Fig. 10. x direction bending moment across the window pane for test 2 at 2.25 ms. The area covered by the plot is shown in Fig. 8. importantly though, it would be beneficial to introduce the ability of the supports to move in the out-of-plane direction. Whilst the effect was relatively limited here, the movement of the supports already introduced some observable errors in the estimates. Also, the set up used in the experiment employed heavy steel sections as supports, thereby limiting their displacements. In practical glazing applications much lighter aluminium sections would be likely to be used, making their deflection much larger. This phenomenon might also become more relevant should the analysis be taken further to the post-crack phase. A further source of errors is the failure of the glass plies. Whilst for ease of comparison and due to uncertainties in the failure time estimations the deflection results were shown for the first 3 ms, it should be remembered that the glass failure is likely to have happened earlier than this point, by more than a millisecond in Test 3. Whilst the cracking event will not immediately affect the deflection shapes, some deviations will be introduced, causing some of the discrepancies.
The failure time estimates are comparable with the experimental results. The estimated failure times of 2.5, 2.25 and 1.75 ms for Tests 1, 2 and 3 respectively are comparable with the experimental times, which were approximately 3, 4 and 2 ms. Whilst the calculated failures took place sooner, especially in the case of Test 2, there was some uncertainty on the times obtained from the experimental data, as they had to be derived from indirect observations. Most importantly, the extra flexibility of the supports in the experiments would have reduced the stresses in the glass, delaying its ultimate failure.
The reactions also showed the same trend as seen in the experimental cases, though the calculated values differed in some cases significantly from the experimental estimates. This was due to both imprecisions in the analytical solutions, which would be amplified during the reaction calculation process, and the high levels of noise in the experimental data used for the reaction measurements. Del Linz et al. [5] stated that the pre-crack reaction results were not likely to be precise. This, together with the simplifications in the supports assumptions, could justify the differences seen in this analysis. In general though, the reactions were shown to be of a similar magnitude to the post-crack reactions calculated from the experiments. This highlighted that these forces need to be considered during the design process of the glazing supports.
The crack spacing results were qualitatively similar to the experimental results. Hooper [21] stated that more dense cracks were observed in bands parallel to the supports, with even denser cracks joining these bands to the corners. A similar situation can be seen in Fig. 12, lending credibility to the analysis.

Conclusion
An analytical method to calculate the deflections and stresses in laminated glass subjected to blast pressure waves was extended by including further terms of the deflection function to produce improved results for larger explosions. The results produced showed that the method was able to predict the deflections magnitude and shape for the panes considered. The results were in    Fig. 8(a) and an image from a typical test (b) [21]. The areas of higher crack density are similar.
all cases representative of the experimental data and in some cases the differences between the simulation and experimental results were minimal. Therefore the results showed the importance of considering the higher order terms in the deflection function when performing analyses for these cases. The results were then used to find the reactions forces acting on the window supports. Both the bending and the shear forces were considered for this. The results showed a pattern similar to that seen in experiments, though with some differences in magnitude. This could be due both to imprecisions in the analytical model and to the high degree of uncertainty in the experimental data used for the measurements.
An estimate of cracking density was also found and compared with experimental records. The patterns found were similar to those described, showing the cracks to be concentrated in similar areas.
The results produced can effectively be used to account more precisely for the pre-crack phase of the blast event on a glass façade. The ability to calculate reaction forces will enable much more accurate designs, increasing their safety and reducing their cost. The analytical approach developed is a very valuable contribution to engineering companies engaged in risk assessment and glass façade design where simpler-to-apply analytically based models are preferable over full-scale FE analysis which can be time consuming for design assessments.