A numerical method for predicting the deformation of crazed laminated windows under blast loading

The design of laminated glazing for blast resistance is significantly complicated by the post-crack behaviour of glass layers. In this research, a novel numerical method based on a semi-analytical energy model is proposed for the post-crack behaviour of crazed panes. To achieve this, the non-homogenous glass cracks patterns observed in literature experimental and analytical work was taken into consideration. It was assumed that, after the glass crazing, further deformations would occur in the cracked edge areas, whilst the central window surface would remain largely undeformed. Therefore, different internal work expressions were formulated for each zone and were then combined in the overall model. The resulting differential equation was then solved numerically. The results obtained were compared with data from four experimental full-scale blast tests for validation. Three of these blast tests (Tests 1 to 3) were presented previously (Hooper et al., 2012) on 1.5 x 1.2 m laminated glazing samples made up with two 3 mm glass layers and a central 1.52 mm PVB membrane, using a 15 and 30 kg charge masses (TNT equivalent) at 13-16 m stand-off. The fourth blast test (Test 4) was conducted on a larger 3.6 x 2.0 m pane of 13.52 mm thickness, using a 100 kg charge mass (TNT equivalent) at a 17 m stand-off. All blast Reaction forces also with experimental estimates. The predictive ability of the proposed accurate to produced rapidly, improving structures resistance to such loadings. thick PVB membrane. No additional treatments besides the normal lamination process were applied to bond the layers of the sample. A white enamel coating was applied to the inside face to facilitate the DIC instrumentation. The glazing pane was mounted in a mild steel frame using a two part silicone sealant joint. This provided a continuous support which had properties close to a pin condition. The test piece was positioned in a concrete cubicle.


Introduction
The blast resistance of glazing is an important consideration when designing against explosions.
Monolithic annealed glass panes produce dangerous shards due to the inherent low fracture toughness of the material. Fragments are propelled both inside and outside the building space and can cause significant injuries and damage. After fracture, residual blast pressures are able to penetrate the building envelope, causing further injuries to occupants and equipment.
Laminated glazing, comprising layers of glass and Polyvinyl Butyral (PVB) membranes, is significantly more resilient to blast loads [1]. After the glass layers craze, the glass fragments remain bonded to the polymer membrane, which need to be retained in the frame. The PVB membrane can then deform significantly, absorbing large amounts of energy and preventing blast pressures from entering the interior.
The behaviour of the crazed pane is complex to model in detail. To gain understanding of the laminated material, several experimental studies have been performed on glazing panels subject to shock loading, either with blast tests [2][3][4][5][6][7][8][9] and in shock tubes [10]. The results of such tests have been used recently by many researchers to validate finite element analyses (FEA) of both impact and blast loading [2,3,[11][12][13][14][15][16][17]. Whilst these models are often able to predict the behaviour of the system, their use requires significant specialist knowledge and computer time.
Analytical solutions instead can produce relatively rapid results, though they cannot account for the same variety of situations as FEA models. A commonly used analysis method is given by single degree of freedom approximations. These can use constants derived either empirically or from first principles to estimate the windows behaviour and therefore produce structural designs [1,18]. This approach is currently used for the design of elements, as it can produce accurate data based on extensive tests databases [1]. Some authors [19][20][21][22] proposed instead more detailed analytical models, as these could be used to predict additional aspects of the glazing behaviour. Wei and Dharani [23] employed Von Karman large deflection theory to simulate the deformation of window panes. The authors used their results to also calculate the probability of window failures and the likely crack locations [24]. In this analysis the authors used a single cosine function to represent the deflected shape of the window. Their work was expanded by Del Linz et al. [25] employing a higher order deflection function and comparing the analytical results with DIC recorded data. The research used three sets of blast test data (Tests 1 to 3) and showed that the analytical technique considered could accurately predict deflections before the glass failure and likely crack concentrations in the failed panels, as well as the reaction forces imposed on the supporting structures. It therefore could provide additional details on the structural behaviour compared with previous approaches, using potentially a smaller number of initial assumptions.
The condition of a glazed pane after it has crazed is more complex than before the glass failure. The crazed pane material is no longer homogenous as its properties are likely to depend on the density and orientation of the cracks. As shown experimentally [4] and through the pre-crack solutions mentioned above [25], in cases where the blast pressures deform the glass impulsively, the crack pattern tends to be denser along bands corresponding to the higher bending stresses in the pane at the points of initial fracture. This influences the deformation which follows, especially in the case of stronger blast excitations. A significant change in the panes curvature occurs at these higher crack density locations throughout the glass deformation history. This effect is indicated by the arrows in Error! Reference source not found. for Test 3, where a clear change in curvature is visible approximately 1/3 of the distance from the edge to the centre of the pane. Error! Reference source not found. shows the locations of the surface cuts plotted in Error! Reference source not found. and in subsequent figures. In all cases the shortest pane side was assumed to be aligned with the x-axis.  Whilst the detailed physical characteristics of the systems are different, the presence of these lines of high curvature highlight a similarity with the plastic large deflection of plates commonly analysed with yield line theory [26]. This analysis is based on equating the external work performed by the loads to the internal energy stored in the system, which is localised at a series of failure yield lines. For laminated glass panes, the energy absorbing capacity due to out-of-plane bending of the cracked glass is close to zero, therefore a different energy absorption method needs to be assumed. As the deflections of the system cannot be assumed to be small, the membrane forces acting on the pane will be significant and represent a possible mechanism for the development of internal energy. Therefore, whilst in other applications the internal work is calculated using the increase of rotation angle at the yield lines and a plastic bending capacity, in this case the elongation of the membranes was considered. As mentioned above, the glass cracks were shown to be concentrated along the edges of the panes in the blasts tests considered in previous studies. Therefore, for this research, it was assumed that the material deformations due to the membrane forces would also be concentrated in these areas, which would be limited by the window supports and the lines of high curvature highlighted in Error! Reference source not found.. The behaviour of such cracked glass was assumed to be similar to that measured by Hooper et al. [4]. A differential equation solution for the systems of interest was derived equating the external energies applied by the blast to the internal membrane energy caused by the deformation of these highly cracked areas, respecting the assumptions of the behaviour described above.
To validate the proposed analytical methods using data comprising different pane dimensions and aspect ratios, data from four blast tests were used. Three of these (Tests 1 to 3) were obtained from Hooper et al. [4] and have been used previously to validate the pre-crack analytical solution [25]. Additionally, the results from a further test (Test 4) conducted by Hooper are presented here and compared with both precrack and post-crack analytical solutions [27]. The glazing pane used for Test 4 was 3.6 m x 2.0 m. The glass make-up was given by two 6 mm glass layers interlayered with a 1.52 mm PVB layer.
The analytical solution for the post-crack behaviour was therefore compared with the experimental data from four tests. An estimate of the reactions based on the model results was also produced for the first three tests so as to compare these results with those produced from experimental data [28]. Due to possible uncertainties in the derivation of the loading function, a small sensitivity study was also conducted to assess the influence of the fitting paramaters.
The aim of this work was to provide a more flexible tool for designers, which could provide additional details, such as the entire deformation history and reactions forces calculated in this work, compared with previously used design methods. Additionally, as per the single degree of freedom approach, calculations would be relatively rapid, therefore offering a useful tool for the practical design of these structures.

Experimental program
The evaluation procedures for all four tests were very similar. Tests 1 to 3 have been described in Hooper et al. [4]. They were conducted on 1.5 x 1.2 m laminated glazing samples made up with two 3 mm glass layers and a central 1.52 mm PVB membrane. The blast loading was a 15 and 30 kg charge masses (TNT equivalent) at 13-16 m stand-off. Test 4, conducted by Hooper [27] on a pane, 3.6 m x 2.0 m, was aimed at measuring the blast response of larger laminated samples. The blast loading was produced by a 100 kg charge mass (TNT equivalent) at a stand-off of 17 m. The glazing pane make-up was two 6 mm layers of annealed glass bonded to a central 1.52 mm thick PVB membrane. No additional treatments besides the normal lamination process were applied to bond the layers of the sample. A white enamel coating was applied to the inside face to facilitate the DIC instrumentation. The glazing pane was mounted in a mild steel frame using a two part silicone sealant joint. This provided a continuous support which had properties close to a pin condition. The test piece was positioned in a concrete cubicle.
During the test, the DIC technique was employed to collect full field 3D deformations and strains of the back face of the pane, as was done in tests 1 to 3 [4]. To achieve this, a random pattern of black dots ("speckles") was applied over the white enamel. Two high speed cameras were employed to capture simultaneous pictures of this pattern from different angles. GOM Aramis DIC software [29] was used to analyse the captured images and produce the strains and deflections in three dimensions over the whole surface. The software was firstly calibrated using the procedure indicated by the manufacturer. With this information, the program could determine the relative positions of the cameras as well as correct for lens distortions in the images. The speckles on the samples surface were then tracked and the deflections and strains determined following the methods described by Schreier et al. [30]. The test arrangement is shown in Error! Reference source not found.. In this case, Photron type SA3 cameras were employed and data were recorded at a frequency of 2000 Hz.
The pressure gauges employed failed. Therefore the reflected peak pressure and impulse were calculated using Air3D [31] and CONWEP [32] at the gauge location.

Analytical method
The results from Tests 1 to 4 were employed to validate analytical solutions for the deflection of the blast loaded glazing panes. The pre-crack solutions were found following the method described by Del Linz et al. [25], though the solution method is summarised below. A post-crack solution was then developed using a technique based on the yield line theory, as the new condition of the system required a different approach compared with the pre-crack situation. Both these systems required several assumptions to be made regarding the loading and the material properties of the glazing panes both before and after the glass crack.

Dynamic loading
It was decided to represent the blast loading using the same approach as earlier authors [23,25,[33][34][35] by employing the equation [36]: where p(t) is the pressure at time t, p0 is the initial peak, td is the positive phase duration and α is an exponential decay coefficient. The equation can represent both the positive and the negative phases of the blast loading. Whilst in the cases considered in this work the positive pressure phase tended to cause the glazing failure, the flexibility of the pressure representation could be useful for more general application.
The shape of the curve is shown in Error! Reference source not found.. As discussed in previous work [25], the reflected pressure was not experimentally recorded during the first three tests considered.
Therefore in these cases the results of computational fluid dynamics analyses were used to fit the parameters of the equation above. The same procedure was also applied to the data of test 4. The CONWEP based tool included in LS-Dyna [37] was used. Whilst it would have been equivalent to use the original curves on which the tool was based, it was decided that using the inbuilt tool would be more convenient in this study. A simple FEA model was created including the glazing pane. This was modelled as continuously supported for simplicity, as the software simulation did not include fluid-structure interactions. The 2D mesh used was composed of 50 mm square elements. The thickness was assumed to be 13.52 mm. Glass material properties were employed. It should be noted that the purpose of the finite element analysis was solely to determine the blast loading pressures. Therefore, the details of the material model and sections would not influence the results. The blast was modelled using the command LOAD_BLAST_ENHANCED and assuming a surface burst. The blast height was assumed to be at the mid height of the panel, and the experimental stand-off of 17 m was used. The charge weight was also assumed to be 100 kg as per the experiment. The equation given above was fitted to the pressure time history using a simple minimisation of square error technique. The final parameters for all the cases are given in table 1. Data regarding the blast tests set ups, including the charge weight used and the stand-off distances, were also included in the same table. The positive impulse applied was calculated integrating the theoretical results. As the lack of direct measurements caused some uncertainty in the pressure estimates, it was decided to also conduct a small sensitivity study with the proposed models varying the impulse, peak pressure and the positive phase durations. The cases considered are listed in table 2. The central deflections were compared to assess the influence of these parameters.

Material properties
For the pre-crack analysis, the same material properties used in the previous study [25] were employed.
The composite Young's modulus and Poisson ratio were found using the rule of admixture: where E is the equivalent Young's modulus of the composite material, Eg and Ep are the Young's moduli of the glass and PVB membrane respectively, and hg and hp are the thicknesses of the glass and PVB.
where υ is the Poisson's ratio of the composite material and υg and υp are the Poisson's ratios of the glass and PVB layers. The total mass per unit area was found through the formula: where μ is the mass per unit area and ρg and ρp are the densities of the glass and PVB layers. The glass and PVB material constants are listed in table 3. As before, a failure stress of 100 MPa was assumed to determine the failure time of the glass [25]. This value is assumed to take into consideration the strain rate effects magnifying the material capacity during dynamic deformations. The value of 80 MPa assumed by previous authors [1,4] was equivalent to the 90% exceedance failure strength at a loading rate of 2 MPa/s (50 MPa) enhanced by strain rate effects.
However, the results of the pre-crack solution [25] indicated that this value under predicted the panes failure times. As the deformation shapes were accurately predicted, it was considered likely that the material used had a capacity higher than previously assumed. A strength of 100 MPa instead produced realistic failure times. Cormie et al. [38] suggested that the dynamic effects on the material can be represented by equation 5: where t1 and σ1 were the original failure time and stress and t2 and σ2 were the new failure time and stress.
k was a constant given as 16 at temperatures below 150º C. As t2 was approximately 2.5 ms in the tests considered, a failure stress of 100 MPa would require a 2 MPa/s capacity of 56 MPa, still below the median failure limit of the data reported by Cormie et al.. Therefore, the limit used was considered realistic in the cases considered in this study.
The post-crack material behaviour of the glazing pane instead showed more complex, inelastic characteristics during the blast tests and in laboratory experiments [4,[39][40][41]. Specifically, it appeared that during laboratory tensile tests a load plateau was reached after small deformations and that subsequent elongations took place at this stress level [4]. A material model showing this behaviour therefore needed to be applied to the present analysis. Galuppi and Royer-Carfagni [42] proposed an analytical method to calculate an equivalent elastic modulus of the cracked section. Alternatively, in previous work [4,28] a Johnson Cook material law was assumed for FEA models and experimental analysis as it could represent both the plateau stress behaviour and the influence of strain rate on the results. Whilst such a material model could be used for this work, it was decided that its mathematical representation would complicate unduly the differential equations to be derived. Additionally, due to the inherent simplifications needed to derive the analytical solution, incorporating the full Johnson Cook law was not considered necessary. However, as it was important to include a material model which could represent the effect of the strain rate, the relationship presented by Hooper [27] was included. In his work, Hooper suggested that the plateau stress levels reached could be related to the strain rate of the tests through a logarithmic function. This was given by Equation 5: where mp, σP, and 0 &are material constants given in table 4. The fit is shown in Error! Reference source not found.. Therefore, in this work, the crazed panes were assumed to behave in a perfectly plastic manner, with the plateau stress level determined with the calculated strain rate.

Pre-crack solution
The analysis of the pre-crack deformations was performed following the method described by Del Linz et al. [25], which employed a solution of the von Karman large deflection equations. These are given by: where F is the Airy's stress function, E is the Young's modulus of the material, w is the out of plane deflection at all points, D is the plate bending stiffness and μ is the mass per unit area of the structure. The loading function described above was assumed for this stage. The letter subscripts after the commas indicate differentiations of the variable in the space coordinate indicated, such that w,xy is equal to 2 dw dxdy .
Pinned boundary conditions were used: where a and b are the dimensions of the panel and the 0,0 point is in the centre of the system.
To solve the equations, the following deflection shape was assumed, cos cos mn mn mx ny wh ab where h is the total thickness of the panes, m and n are integer odd constants and φmn(t) are a series of constants varying in time characterising the amplitude of each mode of vibration at each time point.
This assumption was employed together with Eq. 7 to derive the Airy's stress function (F): Finally, Eq. 8 was changed as shown in Eq. 12 and the previous results were substituted in. The Galerkin method was applied and the system of equations shown in Eq. 13 was obtained. This could then be solved numerically to obtain the constants φmn and therefore the deflections at all points.
These results were then used to calculate the bending and membrane stresses acting at all points of the pane, as done previously. With this data the failure time of the glass was estimated employing the limiting stress described above. Additionally, a crack density across the pane was calculated making use of the strain energy results. The method derived by Wei and Dharani [24] and utilised by Del Linz et al. [25] was applied for this exercise. In this, the energy changes due to the crack formations are used to equilibrate the internal energy at the time of the cracks. The internal energy in the system is given by the components: where U is the total energy, U0 is the strain energy, Ua is the drop in energy due to a crack being formed and Uγ is the surface energy required to create a crack. The total and the strain energy reduce to 0 when a crack is being formed, giving: where: In this, lc is the length of crack and γs is 3.9 J/m 2 [43]. Δx and Δy are the dimension of each facet. Equating the energies lc for each facet can be found. This can then be used as an indication of the likely crack density at the facet.

Post-crack solution
The post-crack solution was derived representing the deformation of the system by assuming the presence of some yield lines which would act as hinges in the panel. All further deformation was assumed to take place between such locations and the supports, where the glazing presented a higher crack density and therefore a lower stiffness. This implied that the areas in between such lines would not undergo further deformation, potentially limiting the accuracy of the simulations. The implications of this assumption are discussed at length in the discussion section using the results of the validation presented below.
A key requirement of the method was to choose the location of the yield lines. These parameters can influence the results significantly. The yield lines were assumed to be parallel to the window edge and to connect the junctions of these lines to the corners, as shown in Error! Reference source not found.. This was done with reference to the results obtained from both experimental and analytical results obtained previously [25]. The critical parameters of the analysis are the dimensions P and Q. They were assumed to be determined by the location of the highest principal tensile stress derived in the elastic plate analysis. However, it was considered that, whilst this point would indicate the position of the finest cracks, a lower crack density could be sufficient to justify the positioning of the yield line. Therefore, it was decided to site the lines further towards the centre, at the point at which the crack concentration in the pane dropped below a predetermined percentage of the maximum. This percentage was determined through trial and error in test 1 and finally it was decided that a value 60% below the peak crack concentration in each pane provided the best results and was kept constant in all cases.
A basic equation for the full internal and external work equilibrium of a plate was given by Jones [44]: where pi is the externally applied load in horizontal direction i ( The main assumption made in the mathematical derivation was that all the deformations were taken into account as if they took place at the yield lines locations. Therefore, all the area integrations could be assumed not to contribute to the internal work. Whilst, as indicated previously, the deformations in the system were assumed to occur in the areas between the yield lines and the supports, it was found to be convenient to calculate the displacements and hence the energy absorbed at the yield lines locations. The internal energy could then be integrated along such lines, simplifying the required algebra and the final differential equation. As a further simplification, the bending moments at the yield lines were assumed to be 0 and therefore were not considered in the final solution. The deflections ui and the in plane external forces were also assumed to be 0. Therefore Eq. 18 reduced to: If only forces perpendicular to the yield lines are assumed, the right hand side of the equation can be proven to be equivalent to multiplying the membrane force by the rate of deformation of the membrane in the force direction. As this is a more convenient expression, it was used in this research: where l is the length of the membrane.
As shown in Error! Reference source not found., the area of the pane could be split into three zones. In the central area "a", it was assumed that the material moved out of plane by a uniform deflection W.
In the areas "b" and "c" the deflection was assumed to be equal to W at the central yield line location and to decrease linearly to 0 at the support. Substituting this into the equation, the external work for areas "b" and "c" were given respectively by: In the equations above L, Q, B and P are defined as shown in Error! Reference source not found.. Therefore, after some simplification, the total external work and kinetic energy for a quarter area is given by: The internal work needed to be calculated along lines 1, 2 and 3. Along line 1, the rate of change in length of the membrane is given by: The material strain and strain rates can be shown to be equal to: where ε1 is the strain of the material between line 1 and the top support. The stress in the membrane is: With some algebraic manipulation, the internal energy is given by: where z is a coordinate along the line length. A similar derivation for line 2 gives: The derivation of the internal work for line 3 follows the same principle, although the algebra is more complex due to the varying deflection along the yield line and its inclined orientation. It can be proven though that the rate of change of length of the membrane between the line and the top support is: 31  These results were used to derive a strain and strain rate level on both sides of the line. Finally, the stress on the two sides of the yield line was found to be equal to: 31 10 Therefore, the internal work calculated along this line was: 31  This equation was then solved using the material properties and loading described above. Matlab [45] was employed to produce a numerical solution. The pre-crack solution at the failure time was used to obtain the initial conditions of the deflection velocity and acceleration of the pane. It should be noted that the equations above were solved equating the change in external work with the changes in the internal and kinetic energies. As the strain energy accumulated before cracking was assumed to be dissipated with the failure of the glass layers, the system energy was not conserved. It was instead assumed that the internal energy would be zero at the beginning of the post crack deformation phase. This is similar to the assumptions made in single degree of freedom systems, where often the resistance function is assumed to decrease to close to zero at the glass failure deflection, dissipating significant amounts of energy [18].
Once the deflections were found, the reaction forces could also be calculated for comparison with the experimental results. To achieve this, the membrane force of each area was multiplied by the sine of the glass inclination angle to find the out of plane components, which were then plotted together with the data obtained in previous research [28] for Tests 1 to 3.

Experimental
The data collected during Test 4 were analysed as described above. The peak pressure estimated with Air3D was equal to 176 kPa, delivering a total impulse of 773 kPa ms. The CONWEP estimate used for the analytical simulation was as shown in table 1. The data showed that the load deformed the window, producing the highest strain concentrations along lines parallel to the support frame, as was seen in previous tests [4]. These concentrations affected the distribution of the fractures of the glass, as can be seen by the higher strains along such a pattern at later times. The glass panes failed approximately 1.5 ms after the blast wave arrival.
The ultimate failure was caused by a failure of the bond between the glass and the silicone. This was facilitated by the enamel coating, which could have weakened the adhesion between the two materials.

Pre-crack solution
The solution detailed above was applied to Tests 1-3 as shown in previous research [25] and additionally to Test 4. The results of this analysis for the latter test are presented in Error! Reference source not found.. These showed that the technique could be applied successfully to glazing panes of different dimensions and aspect ratios to the ones considered previously. The predicted deflections tended to be close to the experimental records at the centre of the window, with a maximum difference of 2 mm. The shape nearer the supports was not approximated as closely as seen in previous tests. The experimental data suggests that significant support deflections, up to 5 mm, took place at this early stage, which would not have been modelled by the analytical method. The calculated crack concentrations are shown in Error! Reference source not found.. The analysis showed that the peak crack concentrations could be found near the corner of the pane. As was done in previous work [25], the crack concentration was expressed as the total crack length for each 10 mm x 10 mm facet.

Post-crack solution
As described above, the pre-crack results were firstly employed to determine the position of the yield lines in the cases considered. A typical result, in this case for Test 3, is shown in Error! Reference source not found.. The differential equation was then solved numerically for each case, employing the last step of the precrack solution to provide the initial conditions for the displacement and velocities of the pane. Typical final results for Tests 1 to 4 are shown in Error! Reference source not found. and Error! Reference source not found., where the deflections along a cut across the centre of the window pane are compared with experimental deflections at the same time.
The central deflections were also compared to estimate whether the overall dynamic behaviour corresponded to the experimental records with regards to its period and maximum deflections. These results are shown in Error! Reference source not found.. The data showed that in most cases the movements calculated with the present method were comparable with the experimentally recorded values. Test 2, 3 and 4 showed the most accurate results, with relative errors generally below 5%. Test 1 instead showed the least accurate deflections, with a maximum relative error of 25%. In general, most discrepancies were due to the analytical solution underestimating the maximum deflections. In Tests 2 and 3, the panes images were showing clear signs of edge failure at 13 ms and 7 ms respectively. Therefore, especially in the case of the latter test, it is possible that some of the discrepancies in the results were due to this. The only case where the deflections were overestimated was test 4. In this the analytical results showed a higher deceleration than the experimental measurements after 10 ms. The periods of the structures were also similar, with a similar deceleration observed in the experimental and analytical data.  The total out of plane reaction forces were calculated for Tests 1 to 3 so as to compare them with the values from previous studies [28]. The values shown include the total reactions along the four edges of each pane. In all cases the calculated values were of similar magnitude to the experimentally estimated data. The results for Test 1 agreed closely with the experimental estimates. Instead, in Tests 2 and 3 the analytical solution overestimated the reaction forces by up to 10.2 kN. A plateau level was present in all cases, showing that the simulated mechanism was realistic, though the forces predicted could be too high, as shown in Error! Reference source not found.. The more accurate prediction for Test 1 also indicated that the quality of this estimate was not affected by the correlation of the calculated and measured central deflections.
The results of the sensitivity study are shown in Error! Reference source not found.. The curves indicated that all the parameters had an important effect on the model behaviour. Reduction in the impulse and in the positive phase duration in cases 1 and 3 caused a marked reduction in the deformations, whilst when these parameters were kept to a similar or higher level the changes in the central deflection seemed to be smaller.

Discussion
The experimental results showed that the larger pane behaved in a similar manner to the smaller windows considered before. The blast impulse and peak pressures were estimated to be relatively high and it is likely that this caused the bands of cracks to be close to the edge of the window and the post-crack deformations to remain concentrated in this narrower band. The elongated aspect ratio of the window did influence the few cracks which appeared in the central area later, which connected the two longer sides.
The deflection at initial cracking could be estimated to be approximately 25 mm, lower than the previously analysed cases. The final deformation at ultimate failure was instead similar at close to 200 mm, though the failure mechanism was different from those previously observed [4].
Whilst it is unfortunate that no direct pressure measurement could be obtained during the test, the CFD and CONWEP analyses had proved to be accurate in previous cases [4]. Therefore it was decided to use numerical techniques as a data source to fit the simplified reflected pressure time history equation for the analytical solutions. Using this input, the pre-crack analysis displayed a good agreement with the experimental DIC data. The central deflections and the deflected shape were realistic, with the latter showing significant deflection and curvature peaks near the edges. Whilst, as mentioned above, the calculated magnitude of the movement was not as accurate away from the central plateau, the strain concentrations also mirrored the DIC results as can be seen comparing Error! Reference source not found. and Error! Reference source not found..
The post-crack analysis produced realistic estimates of the deflections and reactions forces in most of the cases considered. Whilst all four cases showed higher discrepancies with the experimental results than the pre-crack solution, this could be expected since the complexity of the analysis was significantly increased after the glass crazed and a practical analytical model required a larger set of assumptions than the precrack solution. The main assumption was that all further deformation would take place between the yield line locations and the supports. This implied that the central area of the glass would not deform after initial cracking. A second assumption was that the yield lines would not move throughout the analysis time. The results for Test 1 especially showed that these conditions were not always met. The DIC data indicated that the central area increased its curvature as time passed, and the yield line solution underestimated significantly the final deflection, with an error of 23%. The deflected shapes also showed that the position of the higher curvature point (the "yield line") seemed to vary, shifting towards the centre of the window as time passed. This effect was significantly weaker in the cases where the blast loads were higher, such as Test 2, 3 and 4. In these tests, the DIC results indicated that the window deflected impulsively until ultimate failure was reached, with only small amounts of further deformation taking place in the window centre, and the analytical results were much closer to the experimental data.
The assumed positon of the yield lines could also affect the results significantly. As explained above, this parameter was assumed to be related to the peaks in crack concentrations calculated in the pre-crack phase of the loading. However, it proved difficult to derive a physically based criterion to account for this. Instead, the empirical method described was used. Whilst this represents a less general solution, the limit used in this work was fitted to one case and used for the other blast loading and panes dimensions with good outcomes. It is therefore suggested that the proposed parameter is a good approximation for these common loading and structural details combinations.
The small sensitivity study showed that the impulse and the positive load duration were the key parameters of the pressure function. The positive phase duration was similar to the time required to reach the maximum pane deformation and its failure. Therefore, a small reduction in the load time would have caused the negative loading phase to affect the deflections more significantly, causing the marked reduction in the central deflection seen in cases 1 and 3. Instead, whilst the changes in the impulse and pressure in cases 2 and 4 affected the overall results somewhat, these effects were much less significant, especially at later times.
The pre-crack results for Test 4 show that the methods currently cannot account for support movement.
However, this limitation seemed to be less important for the post crack solutions, which were accurate also for cases where more significant support movement was observed due to its relatively small magnitude of generally less than 10 mm.
The reactions estimates produced results similar to the experimental values calculated in previous work [28]. It is likely that the discrepancies observed were due to the assumptions made. The reactions could have been especially affected by the yield lines movement in the tests. This meant that in the experiments the angle of the laminate remained approximately constant after reaching a maximum close to 30° [4], whilst in the simulations the inclination continued to increase, producing higher out of plane reaction forces.
It should be noted that the results presented here were validated using blast tests with scaled distances between 3.6 m/kg 1/3 and 5.3 m/kg 1/3 . At this level of loading, the glass panes generally behaved impulsively at the beginning of the loading, causing the observed pattern of cracks to form. This pattern was then used in this work to calculate the yield line position and to assume a deformation mechanism for the system. Should the loading be significantly lower than the levels considered here, it is possible that a more sinusoidal shaped deformation would be present before the glass failure, invalidating the basic assumptions used for the derivation. Additionally, as mentioned above, lower blast levels would also cause relatively more deformation in the central un-cracked area, as observed in the results for Test 1.
However, such lower levels of loading might also represent a less critical case for the glazing, as its ultimate failure would be less likely. It is therefore suggested that the methods presented here, whilst unable to simulate the behaviour of all loading levels in their present state, could be useful for a critical range of blast pressures which can cause failures in the glazing system. The results indicated that these distances should not exceed significantly 5 m/kg 1/3 . With the available data it is more difficult to establish a lower limit, though it is likely that the system would be accurate at values lower than 3.6 m/kg 1/3 . Further research could be performed to improve on these issues. Ideally, a solution allowing the yield line positions to move should be developed, as this would resolve the issues encountered here. However, the present solution can already provide a powerful and fast tool to analyse laminated glazing panels up to their ultimate failure, providing guidance for the design of structures.

Conclusions
A yield line approximation was derived to model the post-crack behaviour of laminated glazing loaded by blast. Membrane forces were considered dominant for this phase of the deformation and the standard derivations available in literature had to be modified to account for this. Results from a previous analytical solution covering the pre-crack phase of deformations were employed to obtain the starting conditions of the panel after glass failure as well as likely positions of the yield lines. This new method was then applied to four blast tests to validate its results. These included the three blast tests considered when developing the pre-crack analytical solution and a new experiment conducted on a panel with a different aspect ratio and thicker make-up. The aim of this was also to verify whether such panes would behave similarly to the smaller units tested previously. The experimental results suggested that this was the case. The DIC technique was able to capture the deflections and strains, and hence the likely crack locations, well and the data could be used to validate the proposed analytical solutions. The ultimate failure mode of the pane also showed the importance of the bond between the glass and the supports, as in this case the pane capacity is likely to have been limited by the rebate joint failure. This was probably caused by the enamel coating employed, highlighting the effect of relatively minor aspects of glazing design.
The analytical results indicated that the proposed method is able to represent the post-crack deformation of the panels. The central deflections were closely matched in three of the four tests, whilst the reactions magnitudes and time history behaviour were approximately reproduced in all cases. The results did though highlight some limitations in the proposed approach. One of the basic assumptions made was that all the deformation after the glass failure would take place at the yield lines, which was not the case in Test 1. Additionally, the yield lines were assumed not to move, which is likely to have caused the discrepancies in the reaction forces seen in Tests 2 and 3.
There is scope to incorporate the possibility of deformations in the central, non-cracked area and for the yield lines to shift. However, the present results have produced realistic estimates of the system behaviour for this complex physical problem. It is, therefore, hoped that the method will be useful to produce safer and more economical glazing designs under these extreme conditions.